Final version of the bitant-to-fullmode converter. Now it also flips the sign

of reflected data if called with '-flipsign'.
Also added utility program to merge HDF5 data files into a single one
(useful to combine several timesteps or different variables).


git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGHIO/IOHDF5/trunk@90 4825ed28-b72c-4eae-9704-e50c059e567d

2 files changed, 427 insertions, 63 deletions

diff --git a/src/util/hdf5_bitant_to_fullmode.c b/src/util/hdf5_bitant_to_fullmode.c
index 09090d2..bafe9ff 100644
--- a/src/util/hdf5_bitant_to_fullmode.c
+++ b/src/util/hdf5_bitant_to_fullmode.c
@@ -3,7 +3,8 @@
    @date      Thu 10 Jan 2002
    @author    Thomas Radke
    @desc
-              This utility program turns bitant mode datasets into full mode.
+              This utility program copies a Cactus HDF5 datafile turning bitant
+              mode datasets into full mode.
    @enddesc
    @version   $Id$
  @@*/
@@ -11,19 +12,14 @@
 #include "cctk.h"
 
 #include <hdf5.h>
-#include <errno.h> 
 #include <stdio.h> 
 #include <stdlib.h>
 #include <string.h>
-#include <sys/stat.h>
-#ifdef HAVE_UNISTD_H
-#include <unistd.h>     /* sysconf(3) */
-#endif
 #include "BetaThorns/IOHDF5Util/src/ioHDF5UtilGH.h"
 
 /* the rcs ID and its dummy function to use it */
 static const char *rcsid = "$Header$";
-CCTK_FILEVERSION(BetaThorns_IOHDF5_hdf5_bitant_to_fullmode_c)
+CCTK_FILEVERSION(AlphaThorns_IOHDF5_hdf5_bitant_to_fullmode_c)
 
 
 /*****************************************************************************/
@@ -60,6 +56,7 @@ CCTK_FILEVERSION(BetaThorns_IOHDF5_hdf5_bitant_to_fullmode_c)
          a single user-supplied argument */
 static char *pathname = NULL;      /* pathname of the current object */
 static unsigned int nerrors = 0;   /* global error counter */
+static int flipsign = 0;           /* flip sign of reflected elements */
 
 /*****************************************************************************/
 /*                           local function prototypes                       */
@@ -96,69 +93,59 @@ static herr_t CopyAttribute (hid_t src, const char *attr_name, void *arg);
 @@*/
 int main (int argc, char *argv[])
 {
+  int i;
   hid_t infile, outfile;
-#if 0
-  hid_t group, attr, dataspace;
-#endif
 
 
+  /* check for the '-flipsign' flag */
+  flipsign = 0;
+  if (argc > 1 && strcmp (argv[1], "-flipsign") == 0)
+  {
+    flipsign = 1;
+    for (i = 2; i < argc; i++)
+    {
+      argv[i-1] = argv[i];
+    }
+    argc--;
+  }
+
   /* give some help if called with incorrect number of parameters */
   if (argc != 3)
   {
-    fprintf (stderr, "Usage: %s <bitant_infile> <fullmode_outfile>\n", argv[0]);
+    fprintf (stderr, "Usage: %s [-flipsign] <bitant_infile> <fullmode_outfile>"
+                     "\n", argv[0]);
     fprintf (stderr, "   eg, %s alp.h5 alp_fullmode.h5\n\n", argv[0]);
     return (0);
   }
 
-  /* open the bitant mode file */
-  infile = H5Fopen (argv[1], H5F_ACC_RDONLY, H5P_DEFAULT);
-  if (infile < 0)
+  H5E_BEGIN_TRY
   {
-    fprintf (stderr, "ERROR: Cannot open input file '%s' as an HDF5 file !\n\n",
-             argv[1]);
-    return (-1);
-  }
+    /* open the bitant mode file */
+    infile = H5Fopen (argv[1], H5F_ACC_RDONLY, H5P_DEFAULT);
+    if (infile < 0)
+    {
+      fprintf (stderr, "ERROR: Cannot open HDF5 input file '%s' !\n\n",argv[1]);
+      return (-1);
+    }
 
-  /* create output file */
-  outfile = H5Fcreate (argv[2], H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
-  if (outfile < 0)
-  {
-    fprintf (stderr, "ERROR: Cannot create output file '%s' !\n\n", argv[2]);
-    return (-1);
-  }
+    /* create output file */
+    outfile = H5Fcreate (argv[2], H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+    if (outfile < 0)
+    {
+      fprintf (stderr, "ERROR: Cannot create HDF5 output file '%s' !\n\n",
+               argv[2]);
+      return (-1);
+    }
+  } H5E_END_TRY
+
+  printf ("\n  ------------------------------------------\n"
+            "  Cactus 4 HDF5 bitant-to-fullmode Converter\n"
+            "  ------------------------------------------\n");
 
-  /* do the recombination by iterating over all objects
-     in the (first) input file written by processor 0 */
+  /* do the copying by iterating over all objects */
   pathname = "";
   CHECK_ERROR (H5Giterate (infile, "/", NULL, CopyObject, &outfile));
 
-#if 0
-  /* now reset the file info attributes in the GLOBAL_ATTRIBUTES_GROUP group
-     to indicate unchunked file data */
-  CHECK_ERROR (group = H5Gopen (outfile, GLOBAL_ATTRIBUTES_GROUP));
-  CHECK_ERROR (H5Adelete (group, "unchunked"));
-  CHECK_ERROR (H5Adelete (group, "nprocs"));
-  CHECK_ERROR (H5Adelete (group, "ioproc_every"));
-
-  unchunked = 1;
-  nprocs = ioproc_every = 1;
-  CHECK_ERROR (dataspace = H5Screate (H5S_SCALAR));
-  CHECK_ERROR (attr = H5Acreate (group, "unchunked", H5T_NATIVE_INT,
-                                 dataspace, H5P_DEFAULT));
-  CHECK_ERROR (H5Awrite (attr, H5T_NATIVE_INT, &unchunked));
-  CHECK_ERROR (H5Aclose (attr));
-  CHECK_ERROR (attr = H5Acreate (group, "nprocs", H5T_NATIVE_INT,
-                                 dataspace, H5P_DEFAULT));
-  CHECK_ERROR (H5Awrite (attr, H5T_NATIVE_INT, &nprocs));
-  CHECK_ERROR (H5Aclose (attr));
-  CHECK_ERROR (attr = H5Acreate (group, "ioproc_every", H5T_NATIVE_INT,
-                                 dataspace, H5P_DEFAULT));
-  CHECK_ERROR (H5Awrite (attr, H5T_NATIVE_INT, &ioproc_every));
-  CHECK_ERROR (H5Aclose (attr));
-  CHECK_ERROR (H5Sclose (dataspace));
-  CHECK_ERROR (H5Gclose (group));
-#endif
-
   /* finally, close all open files */
   CHECK_ERROR (H5Fclose (infile));
   CHECK_ERROR (H5Fclose (outfile));
@@ -188,7 +175,8 @@ int main (int argc, char *argv[])
    @desc
                Iterator recursively called by H5Giterate() for every object
                in the input file
-               It copies the current object or invokes the recombiner on it.
+               It copies the current object or, if the object is a 3D dataset,
+               reflects the data in the z-plane before copying.
    @enddesc
 
    @calls      CopyAttribute
@@ -222,7 +210,7 @@ static herr_t CopyObject (hid_t from,
   hid_t to, datatype, dataspace_from, dataspace_to;
   H5G_stat_t objectinfo;
   char *current_pathname;
-  int i, offset, xyplane_size, ndims, is_simple, retval;
+  int i, offset, xyplane_size, ndims, is_simple;
   size_t objectsize, datatypesize;
   hsize_t dims[4];
   char *data;
@@ -235,7 +223,6 @@ static herr_t CopyObject (hid_t from,
   sprintf (pathname, "%s/%s", current_pathname, objectname);
 
   /* get the output object identifier */
-  retval = 0;
   to = *(hid_t *) _to;
 
   /* check the type of the current object */
@@ -254,9 +241,6 @@ static herr_t CopyObject (hid_t from,
   }
   else if (objectinfo.type == H5G_DATASET)
   {
-    /* current object is an unchunked dataset - copy it as is */
-    printf ("  copying dataset '%s'\n", pathname);
-
     CHECK_ERROR (from = H5Dopen (from, objectname));
     CHECK_ERROR (datatype = H5Dget_type (from));
     CHECK_ERROR (dataspace_from = H5Dget_space (from));
@@ -264,6 +248,7 @@ static herr_t CopyObject (hid_t from,
     CHECK_ERROR (ndims = H5Sget_simple_extent_ndims (dataspace_from));
     if (is_simple && ndims == 3)
     {
+      printf ("  reflecting dataset '%s'\n", pathname);
       CHECK_ERROR (H5Sget_simple_extent_dims (dataspace_from, dims, NULL));
       /* save the bitant zmax and compute the fullgrid zmax */
       dims[3] = dims[0];
@@ -272,6 +257,7 @@ static herr_t CopyObject (hid_t from,
     }
     else
     {
+      printf ("  copying dataset '%s'\n", pathname);
       CHECK_ERROR (dataspace_to = H5Scopy (dataspace_from));
     }
     CHECK_ERROR (to = H5Dcreate (to, objectname, datatype, dataspace_to,
@@ -300,13 +286,29 @@ static herr_t CopyObject (hid_t from,
       if (is_simple && ndims == 3)
       {
         /* copy the +z xy planes into the -z xy planes */
-        for (i = 1; i < (int) dims[3]; i++)
+        for (i = 2; i < (int) dims[3]; i++)
         {
-          memcpy (data + offset - i*xyplane_size,
+          memcpy (data + offset - (i-1)*xyplane_size,
                   data + offset + i*xyplane_size,
                   xyplane_size);
         }
       }
+      /* flip the sign of the reflected elements */
+      if (flipsign)
+      {
+        if (H5Tequal (datatype, H5T_NATIVE_DOUBLE) > 0)
+        {
+          for (i = offset/sizeof (double) - 1; i >= 0; i--)
+          {
+            ((double *) data)[i] = -((double *) data)[i];
+          }
+        }
+        else
+        {
+          fprintf (stderr, "Warning: can only flip sign for datasets of "
+                           "datatype 'native double'\n");
+        }
+      }
       CHECK_ERROR (H5Dwrite (to, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,data));
       free (data);
     }
@@ -329,7 +331,7 @@ static herr_t CopyObject (hid_t from,
   free (pathname);
   pathname = current_pathname;
 
-  return (retval);
+  return (0);
 }
 
 
diff --git a/src/util/hdf5_merge.c b/src/util/hdf5_merge.c
new file mode 100644
index 0000000..3614bc3
--- /dev/null
+++ b/src/util/hdf5_merge.c
@@ -0,0 +1,362 @@
+ /*@@
+   @file      hdf5_merge.c
+   @date      Thu 10 Jan 2002
+   @author    Thomas Radke
+   @desc
+              This utility program takes a list of Cactus HDF5 datafiles,
+              merges them at the group hierarchy level and dumps the resulting
+              tree to a new HDF5 file.
+   @enddesc
+   @version   $Id$
+ @@*/
+
+#include "cctk.h"
+
+#include <hdf5.h>
+#include <stdio.h> 
+#include <stdlib.h>
+#include <string.h>
+#include "BetaThorns/IOHDF5Util/src/ioHDF5UtilGH.h"
+
+/* the rcs ID and its dummy function to use it */
+static const char *rcsid = "$Header$";
+CCTK_FILEVERSION(AlphaThorns_IOHDF5_merge_c)
+
+
+/*****************************************************************************/
+/*                           macro definitions                               */
+/*****************************************************************************/
+/* uncomment the following to get some debugging output */
+/* #define IOHDF5_DEBUG  1 */
+
+/* macro to do an HDF5 call, check its return code, and print a warning
+   in case of an error */
+#define CHECK_ERROR(hdf5_call)                                                \
+          do                                                                  \
+          {                                                                   \
+            int _error_code = hdf5_call;                                      \
+                                                                              \
+                                                                              \
+            if (_error_code < 0)                                              \
+            {                                                                 \
+              fprintf (stderr, "WARNING: line %d: HDF5 call '%s' returned "   \
+                               "error code %d\n",                             \
+                                __LINE__, #hdf5_call, _error_code);           \
+              nerrors++;                                                      \
+            }                                                                 \
+          } while (0)
+
+
+/*****************************************************************************/
+/*                           global variables                                */
+/*****************************************************************************/
+/* NOTE: although it isn't good programming practice
+         we make these variables global for convenience
+         since they are accessed from recursively or
+         indirectly called routines which only get passed
+         a single user-supplied argument */
+static char *pathname = NULL;      /* pathname of the current object */
+static unsigned int nerrors = 0;   /* global error counter */
+
+/*****************************************************************************/
+/*                           local function prototypes                       */
+/*****************************************************************************/
+static herr_t CopyObject (hid_t copy_from, const char *objectname, void *arg);
+static herr_t CopyAttribute (hid_t src, const char *attr_name, void *arg);
+
+
+ /*@@
+   @routine    main
+   @date       Sat 24 Feb 2001
+   @author     Thomas Radke
+   @desc
+               Main routine of the HDF5 file merger
+   @enddesc
+
+   @calls      CopyObject
+
+   @var        argc
+   @vdesc      number of command line arguments
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        argv
+   @vdesc      command line arguments
+   @vtype      char *[]
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               0 for success, negative return values indicate an error
+   @endreturndesc
+@@*/
+int main (int argc, char *argv[])
+{
+  int i;
+  hid_t *infiles, outfile;
+
+
+  /* give some help if called with incorrect number of parameters */
+  if (argc < 3)
+  {
+    fprintf (stderr, "Usage: %s <infile1> [<infile2> ...] <outfile>\n",argv[0]);
+    fprintf (stderr, "   eg, %s alp.time*.h5 alp.h5"
+                     "\n\n", argv[0]);
+    return (0);
+  }
+
+  H5E_BEGIN_TRY
+  {
+    /* open the input files */
+    infiles = (hid_t *) malloc ((argc - 2) * sizeof (hid_t));
+    for (i = 0; i < argc - 2; i++)
+    {
+      infiles[i] = H5Fopen (argv[i + 1], H5F_ACC_RDONLY, H5P_DEFAULT);
+      if (infiles[i] < 0)
+      {
+        fprintf (stderr, "ERROR: Cannot open HDF5 input file '%s' !\n\n",
+                 argv[i + 1]);
+        return (-1);
+      }
+    }
+
+    /* create output file */
+    outfile = H5Fcreate (argv[argc-1], H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+    if (outfile < 0)
+    {
+      fprintf (stderr, "ERROR: Cannot create HDF5 output file '%s' !\n\n",
+               argv[argc - 1]);
+      return (-1);
+    }
+  } H5E_END_TRY
+
+  printf ("\n  -------------------------\n"
+            "  Cactus 4 HDF5 File Merger\n"
+            "  -------------------------\n");
+
+  /* do the copying by iterating over all objects */
+  for (i = 0; i < argc - 2; i++)
+  {
+    printf ("\n  Merging objects from input file '%s'\n", argv[i + 1]);
+    pathname = "";
+    CHECK_ERROR (H5Giterate (infiles[i], "/", NULL, CopyObject, &outfile));
+  }
+
+  /* finally, close all open files */
+  for (i = 0; i < argc - 2; i++)
+  {
+    CHECK_ERROR (H5Fclose (infiles[i]));
+  }
+  CHECK_ERROR (H5Fclose (outfile));
+  free (infiles);
+
+  /* report status */
+  if (nerrors == 0)
+  {
+    printf ("\n\n   *** All input files successfully merged. ***\n\n");
+  }
+  else
+  {
+    fprintf (stderr, "\n\n   *** WARNING: %d errors occured during "
+                     "file merging. ***\n\n", nerrors);
+  }
+
+  return (0);
+}
+
+
+/*****************************************************************************/
+/*                           local routines                                  */
+/*****************************************************************************/
+ /*@@
+   @routine    CopyObject
+   @date       Sat 24 Feb 2001
+   @author     Thomas Radke
+   @desc
+               Iterator recursively called by H5Giterate() for every object
+               in the input file
+               It copies the current object to the output file if it didn't
+               already exist there.
+   @enddesc
+
+   @calls      CopyAttribute
+
+   @var        from
+   @vdesc      identifier for the group the current object belongs to
+   @vtype      hid_t
+   @vio        in
+   @endvar
+   @var        objectname
+   @vdesc      name of the current object
+   @vtype      const char *
+   @vio        in
+   @endvar
+   @var        _to
+   @vdesc      user-supplied argument indicating the output object identifier
+   @vtype      hid_t
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               0 - continue the iteration for following group objects
+               1 - short-curcuit, no further iteration of this group
+   @endreturndesc
+@@*/
+static herr_t CopyObject (hid_t from,
+                          const char *objectname,
+                          void *_to)
+{
+  int i;
+  hid_t to, datatype, dataspace;
+  H5G_stat_t objectinfo;
+  char *current_pathname;
+  size_t objectsize;
+  void *data;
+
+
+  /* get the output object identifier */
+  to = *(hid_t *) _to;
+
+  /* check whether an object by that name already exists */
+  H5E_BEGIN_TRY
+  {
+    i = H5Gget_objinfo (to, objectname, 0, &objectinfo) >= 0;
+  } H5E_END_TRY
+  if (i)
+  {
+    printf ("   object '%s/%s' will not be copied (already exists)\n",
+            pathname, objectname);
+    return (0);
+  }
+
+  /* build the full pathname for the current to object to process */
+  current_pathname = pathname;
+  pathname = (char *) malloc (strlen (current_pathname) +
+                              strlen (objectname) + 2);
+  sprintf (pathname, "%s/%s", current_pathname, objectname);
+
+  /* check the type of the current object */
+  CHECK_ERROR (H5Gget_objinfo (from, objectname, 0, &objectinfo));
+  if (objectinfo.type == H5G_GROUP)
+  {
+    printf ("   copying group '%s'\n", pathname);
+
+    CHECK_ERROR (from = H5Gopen (from, objectname));
+    CHECK_ERROR (to = H5Gcreate (to, objectname, 0));
+    /* iterate over all objects in the (first) input file */
+    CHECK_ERROR (H5Giterate (from, ".", NULL, CopyObject, &to));
+    CHECK_ERROR (H5Aiterate (from, NULL, CopyAttribute, &to));
+    CHECK_ERROR (H5Gclose (to));
+    CHECK_ERROR (H5Gclose (from));
+  }
+  else if (objectinfo.type == H5G_DATASET)
+  {
+    printf ("   copying dataset '%s'\n", pathname);
+
+    CHECK_ERROR (from = H5Dopen (from, objectname));
+    CHECK_ERROR (datatype = H5Dget_type (from));
+    CHECK_ERROR (dataspace = H5Dget_space (from));
+    CHECK_ERROR (to = H5Dcreate (to, objectname, datatype, dataspace,
+                                 H5P_DEFAULT));
+    objectsize = H5Dget_storage_size (from);
+    if (objectsize > 0)
+    {
+      data = malloc (objectsize);
+      CHECK_ERROR (H5Dread (from, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                            data));
+      CHECK_ERROR (H5Dwrite (to, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,data));
+      free (data);
+    }
+    CHECK_ERROR (H5Aiterate (from, NULL, CopyAttribute, &to));
+    CHECK_ERROR (H5Dclose (to));
+    CHECK_ERROR (H5Dclose (from));
+    CHECK_ERROR (H5Sclose (dataspace));
+    CHECK_ERROR (H5Tclose (datatype));
+  }
+  else
+  {
+    fprintf (stderr, "WARNING: Found object '%s' which is not neither a "
+                     "group nor a dataset ! Object will not be copied.\n",
+                     pathname);
+    nerrors++;
+  }
+
+  /* reset the pathname */
+  free (pathname);
+  pathname = current_pathname;
+
+  return (0);
+}
+
+
+ /*@@
+   @routine    CopyAttribute
+   @date       Sat 24 Feb 2001
+   @author     Thomas Radke
+   @desc
+               Iterator recursively called by H5Aiterate() for every attribute
+               of an object (dataset or group)
+   @enddesc
+
+   @var        from
+   @vdesc      identifier for the group or dataset to read the attribute from
+   @vtype      hid_t
+   @vio        in
+   @endvar
+   @var        attrname
+   @vdesc      name of the current attribute
+   @vtype      const char *
+   @vio        in
+   @endvar
+   @var        _to
+   @vdesc      user-supplied argument indicating the group or dataset
+               to copy the attribute to
+   @vtype      hid_t
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               0 - continue the iteration for following attributes
+   @endreturndesc
+@@*/
+static herr_t CopyAttribute (hid_t from,
+                             const char *attrname,
+                             void *_to)
+{
+  hid_t attr, datatype, dataspace, to;
+  size_t attrsize;
+  void *value;
+
+
+  /* get the target group/dataset identifier */
+  to = *(hid_t *) _to;
+
+  /* open the attribute given by its name, get type, dataspace, and value
+     and just copy it */
+  CHECK_ERROR (attr = H5Aopen_name (from, attrname));
+  CHECK_ERROR (datatype = H5Aget_type (attr));
+  CHECK_ERROR (dataspace = H5Aget_space (attr));
+  attrsize = H5Tget_size (datatype);
+  if (H5Sis_simple (dataspace) > 0)
+  {
+    attrsize *= H5Sget_simple_extent_npoints (dataspace);
+  }
+  if (attrsize > 0)
+  {
+    value = malloc (attrsize);
+    CHECK_ERROR (H5Aread (attr, datatype, value));
+    CHECK_ERROR (H5Aclose (attr));
+    CHECK_ERROR (attr = H5Acreate (to, attrname, datatype, dataspace,
+                                    H5P_DEFAULT));
+    CHECK_ERROR (H5Awrite (attr, datatype, value));
+    free (value);
+  }
+  CHECK_ERROR (H5Aclose (attr));
+  CHECK_ERROR (H5Sclose (dataspace));
+  CHECK_ERROR (H5Tclose (datatype));
+
+  return (0);
+}