path: root/src/util/hdf5_bitant_to_fullmode.c

 /*@@
   @file      hdf5_bitant_to_fullmode.c
   @date      Thu 10 Jan 2002
   @author    Thomas Radke
   @desc
              This utility program copies a Cactus HDF5 datafile turning bitant
              mode datasets into full mode.
   @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_hdf5_bitant_to_fullmode_c)


/*****************************************************************************/
/*                           macro definitions                               */
/*****************************************************************************/
/* uncomment the following to get some debugging output */
/* #define IOHDF5_DEBUG  1 */

/* default stencil width if no '-stencil' option was given */
#define STENCIL          2

/* 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 */
static int stencil = STENCIL;      /* stencil width for reflection */

/*****************************************************************************/
/*                           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);
static void FixBoundingBox (hid_t object);


 /*@@
   @routine    main
   @date       Thu 10 Jan 2002
   @author     Thomas Radke
   @desc
               Main routine of the HDF5 bitant-to-full converter
   @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 infile, outfile;


  /* check for the '-stencil' option */
  if (argc > 1 && strncmp (argv[1], "-stencil=", 9) == 0)
  {
    stencil = atoi (argv[1] + 9);
    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 [-stencil=<N>] <bitant_infile> "
                     "<fullmode_outfile>\n", argv[0]);
    fprintf (stderr, "   eg, %s alp.h5 alp_fullmode.h5\n\n", argv[0]);
    return (0);
  }

  H5E_BEGIN_TRY
  {
    /* 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 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"
            "  Using stencil width of %d for reflection...\n", stencil);

  /* do the copying by iterating over all objects */
  pathname = "";
  CHECK_ERROR (H5Giterate (infile, "/", NULL, CopyObject, &outfile));

  /* finally, close all open files */
  CHECK_ERROR (H5Fclose (infile));
  CHECK_ERROR (H5Fclose (outfile));

  /* report status */
  if (nerrors == 0)
  {
    printf ("\n\n   *** All Cactus data successfully converted. ***\n\n");
  }
  else
  {
    fprintf (stderr, "\n\n   *** WARNING: %d errors occured during "
                     "conversion. ***\n\n", nerrors);
  }

  return (0);
}


/*****************************************************************************/
/*                           local routines                                  */
/*****************************************************************************/
 /*@@
   @routine    CopyObject
   @date       Thu 10 Jan 2002
   @author     Thomas Radke
   @desc
               Iterator recursively called by H5Giterate() for every object
               in the input file
               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

   @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)
{
  hid_t to, datatype, dataspace_from, dataspace_to;
  H5G_stat_t objectinfo;
  char *current_pathname;
  int i, offset, xyplane_size, ndims, is_simple;
  size_t objectsize, datatypesize;
  hsize_t dims[4];
  char *data;


  /* 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);

  /* get the output object identifier */
  to = *(hid_t *) _to;

  /* 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));
    FixBoundingBox (to);
    CHECK_ERROR (H5Gclose (to));
    CHECK_ERROR (H5Gclose (from));
  }
  else if (objectinfo.type == H5G_DATASET)
  {
    CHECK_ERROR (from = H5Dopen (from, objectname));
    CHECK_ERROR (datatype = H5Dget_type (from));
    CHECK_ERROR (dataspace_from = H5Dget_space (from));
    CHECK_ERROR (is_simple = H5Sis_simple (dataspace_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] - stencil;
      dims[0] = 2 * dims[3];
      CHECK_ERROR (dataspace_to = H5Screate_simple (3, dims, NULL));
    }
    else
    {
      printf ("  copying dataset '%s'\n", pathname);
      CHECK_ERROR (dataspace_to = H5Scopy (dataspace_from));
    }
    CHECK_ERROR (to = H5Dcreate (to, objectname, datatype, dataspace_to,
                                 H5P_DEFAULT));
    datatypesize = objectsize = H5Tget_size (datatype);
    xyplane_size = 0;
    if (is_simple)
    {
      objectsize *= H5Sget_simple_extent_npoints (dataspace_from);
    }
    if (objectsize > 0)
    {
      if (is_simple && ndims == 3)
      {
        xyplane_size = dims[1] * dims[2] * datatypesize;
        offset = (dims[3] - stencil) * xyplane_size;
        objectsize += offset;
      }
      else
      {
        offset = 0;
      }
      data = (char *) malloc (objectsize);
      CHECK_ERROR (H5Dread (from, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                            data + offset));
      if (is_simple && ndims == 3)
      {
        /* copy the +z xy planes into the -z xy planes */
        for (i = stencil + 1; i <= (int) dims[3]; i++)
        {
          memcpy (data + offset - (i-stencil)*xyplane_size,
                  data + offset + (i+stencil-1)*xyplane_size,
                  xyplane_size);
        }
      }
#if 0
      /* 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");
        }
      }
#endif
      CHECK_ERROR (H5Dwrite (to, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,data));
      free (data);
    }
    CHECK_ERROR (H5Aiterate (from, NULL, CopyAttribute, &to));
    if (is_simple && ndims == 3)
    {
      FixBoundingBox (to);
    }
    CHECK_ERROR (H5Dclose (to));
    CHECK_ERROR (H5Dclose (from));
    CHECK_ERROR (H5Sclose (dataspace_to));
    CHECK_ERROR (H5Sclose (dataspace_from));
    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       Thu 10 Jan 2002
   @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);
}


 /*@@
   @routine    FixBoundingBox
   @date       Mon 18 Feb 2002
   @author     Thomas Radke
   @desc
               Fix the "global_size", "min_ext", and "origin" attributes
               if found in the object (dataset or group) given.
   @enddesc

   @var        object
   @vdesc      identifier for the group or dataset to fix the attributes
   @vtype      hid_t
   @vio        inout
   @endvar
@@*/
static void FixBoundingBox (hid_t object)
{
  hid_t attr, datatype, dataspace;
  double delta[3], min_ext[3];
  int global_size[3], do_it;


  /* try to read the "delta" and "global_size" attributes */
  H5E_BEGIN_TRY
  {
    attr = H5Aopen_name (object, "delta");
    do_it = attr >= 0;
    if (do_it)
    {
      CHECK_ERROR (datatype = H5Aget_type (attr));
      CHECK_ERROR (dataspace = H5Aget_space (attr));
      if (H5Tget_class (datatype) == H5T_FLOAT &&
          H5Sis_simple (dataspace) > 0 &&
          H5Sget_simple_extent_npoints (dataspace) == 3)
      {
        CHECK_ERROR (H5Aread (attr, H5T_NATIVE_DOUBLE, delta));
      }
      else
      {
        do_it = 0;
      }
      CHECK_ERROR (H5Aclose (attr));
      CHECK_ERROR (H5Sclose (dataspace));
      CHECK_ERROR (H5Tclose (datatype));
    }

    if (do_it)
    {
      attr = H5Aopen_name (object, "global_size");
      do_it = attr >= 0;
    }
    if (do_it)
    {
      CHECK_ERROR (datatype = H5Aget_type (attr));
      CHECK_ERROR (dataspace = H5Aget_space (attr));
      if (H5Tget_class (datatype) == H5T_INTEGER &&
          H5Sis_simple (dataspace) > 0 &&
          H5Sget_simple_extent_npoints (dataspace) == 3)
      {
        CHECK_ERROR (H5Aread (attr, H5T_NATIVE_INT, global_size));
        global_size[2] += global_size[2] - 1;
        CHECK_ERROR (H5Awrite (attr, H5T_NATIVE_INT, global_size));
      }
      else
      {
        do_it = 0;
      }
      CHECK_ERROR (H5Aclose (attr));
      CHECK_ERROR (H5Sclose (dataspace));
      CHECK_ERROR (H5Tclose (datatype));
    }
  } H5E_END_TRY;

  /* now fix the "min_ext" and "origin" attributes */
  if (do_it)
  {
    CHECK_ERROR (attr = H5Aopen_name (object, "min_ext"));
    if (attr >= 0)
    {
      CHECK_ERROR (datatype = H5Aget_type (attr));
      CHECK_ERROR (dataspace = H5Aget_space (attr));
      if (H5Tget_class (datatype) == H5T_FLOAT &&
          H5Sis_simple (dataspace) > 0 &&
          H5Sget_simple_extent_npoints (dataspace) == 3)
      {
        CHECK_ERROR (H5Aread (attr, H5T_NATIVE_DOUBLE, min_ext));
        min_ext[2] -= delta[2] * ((global_size[2] + 1) / 2);
        CHECK_ERROR (H5Awrite (attr, H5T_NATIVE_DOUBLE, min_ext));
      }
      CHECK_ERROR (H5Aclose (attr));
      CHECK_ERROR (H5Sclose (dataspace));
      CHECK_ERROR (H5Tclose (datatype));
    }

    CHECK_ERROR (attr = H5Aopen_name (object, "origin"));
    if (attr >= 0)
    {
      CHECK_ERROR (datatype = H5Aget_type (attr));
      CHECK_ERROR (dataspace = H5Aget_space (attr));
      if (H5Tget_class (datatype) == H5T_FLOAT &&
          H5Sis_simple (dataspace) > 0 &&
          H5Sget_simple_extent_npoints (dataspace) == 3)
      {
        CHECK_ERROR (H5Aread (attr, H5T_NATIVE_DOUBLE, min_ext));
        min_ext[2] -= delta[2] * ((global_size[2] + 1) / 2);
        CHECK_ERROR (H5Awrite (attr, H5T_NATIVE_DOUBLE, min_ext));
      }
      CHECK_ERROR (H5Aclose (attr));
      CHECK_ERROR (H5Sclose (dataspace));
      CHECK_ERROR (H5Tclose (datatype));
    }
  }
}