path: root/src/RestoreFile.c

 /*@@
   @file      RestoreFile.c
   @date      Thu Jun 18 16:34:59 1998
   @author    Tom Goodale
   @desc 
   Routines to restore variables from a given file
   @enddesc 
   @history
   @hauthor Gabrielle Allen @hdate 19 Oct 1998      
   @hdesc Changed names ready for thorn_IO
   @endhistory
   @version $Id$
 @@*/


#include <stdio.h>
#include <stdlib.h>

#include "cctk.h"
#include "cctk_Parameters.h"
#include "CactusBase/IOUtil/src/ioGH.h"
#include "CactusPUGH/PUGH/src/include/pugh.h"
#include "ioHDF5GH.h"


/* MPI tag base */
#define STARTUPBASE 1001        


typedef struct {
  cGH *GH;
  int ioproc;
  int ioproc_every;
  int unchunked;
} iterate_info_t;

typedef struct {
  iterate_info_t *it_info;
  int element_size;
  int iohdf5_type;
#ifdef CCTK_MPI
  MPI_Datatype mpi_type;
#endif
} recover_info_t;


/* local function prototypes */
static herr_t processDataset (hid_t group, const char *datasetname, void *arg);


 /*@@
   @routine    IOHDF5i_RecoverVariables
   @date       Fri Jun 19 09:19:48 1998
   @author     Tom Goodale
   @desc 
   Reads in data from an open HDF5 file.
   Code was originally in StartupReader.
   @enddesc 
   @calledby   IOHDF5_Recover
   @history 
   @endhistory 
   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      cGH
   @vio        in
   @endvar
   @var        file
   @vdesc      HDF5 file info structure
   @vtype      fileinfo_t *
   @vio        in
   @endvar
@@*/
int IOHDF5i_RecoverVariables (cGH *GH, fileinfo_t *file)
{
  DECLARE_CCTK_PARAMETERS
#ifdef CCTK_MPI
  pGH *pughGH;               /* PUGH extension handle */
  int proc;                  /* looper */
  CCTK_INT var_info [3];     /* communication buffer for MPI */
#endif
  iterate_info_t info;       /* info to pass down to the iterator routine */


  info.GH = GH;
  info.ioproc = file->ioproc;
  info.unchunked = file->unchunked;
  info.ioproc_every = file->ioproc_every;

#ifdef CCTK_MPI
  /* Get the handle for PUGH extensions */
  pughGH = PUGH_pGH (GH);

  /* Now process the datasets.
     All IO processors read the datasets from their checkpoint file,
     verify their contents and communicate them to the non-IO processors. */

  /* At first the code for the IO processors.
     This holds also for the single processor case. */
  if (CCTK_MyProc (GH) == file->ioproc) {
#endif

    /* iterate over all datasets starting from "/" in the HDF5 file */
    CACTUS_IOHDF5_ERROR (H5Giterate (file->fid, "/", NULL, processDataset,
                         &info));

#ifdef CCTK_MPI
    /* To signal completion to the non-IO processors
       an invalid variable index is broadcasted. */
    var_info [0] = -1;
    for (proc = 1; proc < file->ioproc_every; proc++)
    for (proc = file->ioproc + 1;
         proc < file->ioproc + file->ioproc_every && proc < CCTK_nProcs (GH);
         proc++)
      CACTUS_MPI_ERROR (MPI_Send (var_info, 3, PUGH_MPI_INT, proc,
                        STARTUPBASE, pughGH->PUGH_COMM_WORLD));

  } else {

    /* And here the code for non-IO processors: */
    MPI_Datatype mpi_type;
    MPI_Status ms;
    int index, timelevel, npoints;

    /* They don't know how many datasets there are, because the IO processors
       could skip some on the fly during their consistency checks.
       The IO Processor sends the index of the variable to be processed next.
       So, all non-IO processors execute a loop where the termination condition
       is when an invalid index was received.
    */
    while (1) {
      /* receive the next variable index from my IO processor */
      CACTUS_MPI_ERROR (MPI_Recv (var_info, 3, PUGH_MPI_INT, file->ioproc,
                        STARTUPBASE, pughGH->PUGH_COMM_WORLD, &ms));
      index = var_info [0]; timelevel = var_info [1]; npoints = var_info [2];

      /* check for termination condition */
      if (index < 0)
        break;

      switch (CCTK_VarTypeI (index)) {
        case CCTK_VARIABLE_CHAR:    mpi_type = PUGH_MPI_CHAR; break;
        case CCTK_VARIABLE_INT:     mpi_type = PUGH_MPI_INT;  break;
        case CCTK_VARIABLE_REAL:    mpi_type = PUGH_MPI_REAL; break;
        case CCTK_VARIABLE_COMPLEX: mpi_type = pughGH->PUGH_mpi_complex; break;
        default:
          CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                      "Unsupported datatype %d", CCTK_VarTypeI (index));
          continue;
      }

      /* receive following data from my IO processor */
      CACTUS_MPI_ERROR (MPI_Recv (GH->data [index][timelevel], npoints,mpi_type,
                        file->ioproc, STARTUPBASE,pughGH->PUGH_COMM_WORLD,&ms));
    }
  }
#endif /* MPI */

  return (0);
}


/**************************** local routines ****************************/

/* local routine GetCommonAttributes() reads in the next dataset's attributes
   and verifies them:

  * checks if there is a variable with the name given by the name attribute
  * verifies that this variable still belongs to the same group
  * checks the group data info:
    - group type
    - variable type
    - ntimelevels
    - sizes (rank, dimensions) according to chunking mode

 If there is a mismatch a warning (warning level 2) is printed and value of 1
 is returned to indicate that this dataset should be ignored.
 If successful, the global variable index, the group type and the timelevel
 to restore are stored in {*index, *gtype, *timelevel}, and 0 is returned.
*/

static int GetCommonAttributes (cGH *GH, hid_t dataset, const char *datasetname,
                                int unchunked, int *index, int *grouptype,
                                int *timelevel, int is_group)
{
  pGExtras *extras;
  cGroup groupdata;
  int i, flag, *dims;
  hid_t datatype, dataspace;
  hsize_t rank_stored, *dims_stored;
  int grouptype_stored, ndims_stored, numtimelevels_stored;
  char *groupname, groupname_stored [128], fullvarname [128];

  /* decompose the datasetname, ignore the iteration number */
  if (sscanf (datasetname, "%[^@]@%d@%d", fullvarname, &i, timelevel) != 3) {
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Cannot parse datasetname '%s'", datasetname);
    return (1);
  }

  /* check if there is a matching variable */
  *index = CCTK_VarIndex (fullvarname);
  if (*index < 0) {
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "No matching variable found for '%s'", fullvarname);
    return (1);
  }

  /* read and verify the group name */
  READ_ATTRIBUTE (dataset, "groupname", H5T_C_S1, groupname_stored);
  groupname = CCTK_GroupNameFromVarI (*index);
  if (! CCTK_Equals (groupname_stored, groupname)) {
    CCTK_WARN (2, "Groupnames don't match");
    return (1);
  }
  free (groupname);

  /* verify group type, variable type, dims, sizes and ntimelevels */
  READ_ATTRIBUTE (dataset, "grouptype", H5T_NATIVE_INT, &grouptype_stored);
  READ_ATTRIBUTE (dataset, "ntimelevels", H5T_NATIVE_INT,&numtimelevels_stored);

  if (CCTK_GroupData (CCTK_GroupIndex (groupname_stored), &groupdata) != 0) {
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Cannot get group data for '%s'", fullvarname);
    return (1);
  }
  if (groupdata.grouptype != grouptype_stored) {
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Group types don't match for '%s'", fullvarname);
    return (1);
  }
  if (groupdata.numtimelevels != numtimelevels_stored) {
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Number of timelevels don't match for '%s'", fullvarname);
    return (1);
  }
  /* open the first chunk to determine datatype, dims and sizes
     if the dataset is a chunk group */
  if (is_group)
    CACTUS_IOHDF5_ERROR (dataset = H5Dopen (dataset, "chunk0"));
  CACTUS_IOHDF5_ERROR (datatype = H5Dget_type (dataset));

  /* The CCTK variable type defines do not correlate with the HDF5 defines
     so compare them explicitely here. */
  flag = (H5Tget_class (datatype) == H5T_FLOAT &&
          groupdata.vartype == CCTK_VARIABLE_REAL) ||
         (H5Tget_class (datatype) == H5T_INTEGER &&
          groupdata.vartype == CCTK_VARIABLE_INT) ||
         (H5Tget_class (datatype) == H5T_INTEGER &&
          groupdata.vartype == CCTK_VARIABLE_CHAR) ||
         (H5Tget_class (datatype) == H5T_COMPOUND &&
          groupdata.vartype == CCTK_VARIABLE_COMPLEX);

  CACTUS_IOHDF5_ERROR (H5Tclose (datatype));
  if (! flag) {
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Variable types don't match for '%s'", fullvarname);
    return (1);
  }

  /* verify the dims and sizes */
  CACTUS_IOHDF5_ERROR (dataspace = H5Dget_space (dataset));
  CACTUS_IOHDF5_ERROR (ndims_stored = H5Sget_simple_extent_ndims (dataspace));
  dims_stored = (hsize_t *) malloc (ndims_stored * sizeof (hsize_t));
  CACTUS_IOHDF5_ERROR (rank_stored = H5Sget_simple_extent_dims (dataspace,
                       dims_stored, NULL));
  /* scalars are stored as rank 0 in HDF5 but have rank 1 in Cactus */
  if (rank_stored == 0)
    rank_stored = 1;
  CACTUS_IOHDF5_ERROR (H5Sclose (dataspace));

  flag = 0;
  if (groupdata.dim != rank_stored)
    flag = 1;
  switch (groupdata.grouptype) {
    case CCTK_SCALAR:
      break;
    case CCTK_ARRAY:
    case CCTK_GF:
      extras = ((pGA ***) PUGH_pGH (GH)->variables)[*index][*timelevel]->extras;
      dims = unchunked ? extras->nsize : extras->lnsize;
      for (i = 0; i < groupdata.dim; i++)
        if (dims [groupdata.dim - i - 1] != dims_stored [i])
          flag = 1;
      break;
  }

  free (dims_stored);
  if (flag) {
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Variable sizes don't match for '%s'", fullvarname);
    return (1);
  }

  if (! CCTK_QueryGroupStorageI (GH, CCTK_GroupIndexFromVarI (*index))) {
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Can't read into variable '%s': no storage", fullvarname);
    return (1);
  }

  /* close the first chunk if the dataset is a chunk group */
  if (is_group)
    CACTUS_IOHDF5_ERROR (H5Dclose (dataset));

  *grouptype = groupdata.grouptype;

  return (0);
}


static int IOHDF5_RestoreGS (hid_t dataset, int index, int timelevel,
                             recover_info_t *rec_info)
{
  void *data;
#ifdef CCTK_MPI
  int proc;
  CCTK_INT var_info [3];
  pGH *pughGH;
#endif


  data = CCTK_VarDataPtrI (rec_info->it_info->GH, timelevel, index);

  /* read the data into the local variable ... */
  CACTUS_IOHDF5_ERROR (H5Dread (dataset, rec_info->iohdf5_type, H5S_ALL,
                       H5S_ALL, H5P_DEFAULT, data));

#ifdef CCTK_MPI
  /* ... and communicate it for the MPI case */

  /* Get the handle for PUGH extensions */
  pughGH = PUGH_pGH (rec_info->it_info->GH);

  /* set the variable's index and the timelevel */
  var_info [0] = index; var_info [1] = timelevel; var_info [2] = 1;

  /* send info and data to the non-IO processors */
  for (proc = rec_info->it_info->ioproc + 1;
       proc < rec_info->it_info->ioproc + rec_info->it_info->ioproc_every &&
       proc < CCTK_nProcs (rec_info->it_info->GH);
       proc++) {
    CACTUS_MPI_ERROR (MPI_Send (var_info, 3, PUGH_MPI_INT, proc,
                      STARTUPBASE, pughGH->PUGH_COMM_WORLD));
    CACTUS_MPI_ERROR (MPI_Send (data, 1, rec_info->mpi_type, proc,
                      STARTUPBASE, pughGH->PUGH_COMM_WORLD));
  }
#endif

  return (0);
}


static int IOHDF5_RestoreGA (hid_t dataset, int index, int timelevel,
                             recover_info_t *rec_info)
{
#ifdef CCTK_MPI
  int proc, npoints;
  CCTK_INT var_info [3];
  pGH *pughGH;
  void *buffer, *data;
  hid_t filespace, memspace;
  pGExtras *extras;
#endif


  /* single processor case is easy: just read the whole dataset */
  if (CCTK_nProcs (rec_info->it_info->GH) == 1) {
    CACTUS_IOHDF5_ERROR (H5Dread (dataset, rec_info->iohdf5_type, H5S_ALL,
         H5S_ALL, H5P_DEFAULT, rec_info->it_info->GH->data [index][timelevel]));
    return (0);
  }

#ifdef CCTK_MPI

  /* Get the handle for PUGH extensions */
  pughGH = PUGH_pGH (rec_info->it_info->GH);

  /* Get the pGExtras pointer as a shortcut */
  extras = ((pGA ***) pughGH->variables) [index][timelevel]->extras;

  /* allocate memory for the biggest chunk */
  npoints = extras->rnpoints [rec_info->it_info->ioproc + 1];
  for (proc = 2; proc < rec_info->it_info->ioproc_every; proc++)
    if (npoints < extras->rnpoints [rec_info->it_info->ioproc + proc])
      npoints = extras->rnpoints [rec_info->it_info->ioproc + proc];
  buffer = malloc (npoints * rec_info->element_size);

  /* set the variable's index and timelevel to restore */
  var_info [0] = index; var_info [1] = timelevel;

  /* now loop over the group of processors associated to each IO processor */
  for (proc = rec_info->it_info->ioproc;
       proc < rec_info->it_info->ioproc + rec_info->it_info->ioproc_every &&
       proc < CCTK_nProcs (rec_info->it_info->GH);
       proc++) {

    /* read own data directly into variable */
    if (proc == rec_info->it_info->ioproc)
      data = rec_info->it_info->GH->data [index][timelevel];
    else
      data = buffer;

    if (! rec_info->it_info->unchunked) {
      hid_t chunk;
      char chunkname [32];

      /* Chunked data is stored as separate chunk datasets within a group.
         So open, read and close the separate chunks here. */
      sprintf (chunkname, "chunk%d", proc - rec_info->it_info->ioproc);
      CACTUS_IOHDF5_ERROR (chunk = H5Dopen (dataset, chunkname));
      CACTUS_IOHDF5_ERROR (H5Dread (chunk, rec_info->iohdf5_type, H5S_ALL,
                           H5S_ALL, H5P_DEFAULT, data));
      CACTUS_IOHDF5_ERROR (H5Dclose (chunk));

    } else {
      int i, dim;
      hssize_t *chunk_origin;
      hsize_t *chunk_dims;

      /* get the dimension of the variable */
      dim = CCTK_GroupDimI (CCTK_GroupIndexFromVarI (index));
      chunk_origin = (hssize_t *) malloc (dim * sizeof (hssize_t));
      chunk_dims = (hsize_t *) malloc (dim * sizeof (hsize_t));

      /* Unchunked data is read as one hyperslab per processor.
         So prepare the memspace and the filespace and read the hyperslab. */
      for (i = 0; i < dim; i++) {
        chunk_dims [dim - 1 - i] = extras->rnsize [proc][i];
        chunk_origin [dim - 1 - i] = extras->lb [proc][i];
      }

      CACTUS_IOHDF5_ERROR (filespace = H5Dget_space (dataset));
      CACTUS_IOHDF5_ERROR (memspace = H5Screate_simple (dim, chunk_dims, NULL));
      CACTUS_IOHDF5_ERROR (H5Sselect_hyperslab (filespace, H5S_SELECT_SET,
                           chunk_origin, NULL, chunk_dims, NULL));

      CACTUS_IOHDF5_ERROR (H5Dread (dataset, rec_info->iohdf5_type, memspace,
                           filespace, H5P_DEFAULT, data));

      CACTUS_IOHDF5_ERROR (H5Sclose (memspace));
      CACTUS_IOHDF5_ERROR (H5Sclose (filespace));

      free (chunk_dims);
      free (chunk_origin);
    }

    /* send the index and the data to the non-IO processors */
    if (proc != rec_info->it_info->ioproc) {
      var_info [2] = extras->rnpoints [proc];
      CACTUS_MPI_ERROR (MPI_Send (var_info, 3, PUGH_MPI_INT, proc,
                        STARTUPBASE, pughGH->PUGH_COMM_WORLD));
      CACTUS_MPI_ERROR (MPI_Send (data, extras->rnpoints [proc],
               rec_info->mpi_type, proc, STARTUPBASE, pughGH->PUGH_COMM_WORLD));
    }
  }

  free (buffer);
#endif

  return (0);
}


static herr_t processDataset (hid_t group, const char *datasetname, void *arg)
{
  DECLARE_CCTK_PARAMETERS
  pGH *pughGH;
  ioGH *ioUtilGH;
  ioHDF5GH *myGH;
  int index, gtype, timelevel;
  iterate_info_t *it_info = (iterate_info_t *) arg;
  recover_info_t rec_info;
  hid_t dataset;
  H5G_stat_t object_info;
  int is_group;

  /* skip the global attributes and GH extensions groups */
  if (! strcmp (datasetname, GLOBAL_PARAMETERS_GROUP) ||
      ! strcmp (datasetname, GHEXTENSIONS_GROUP))
    return (0);

  /* Get the handle for PUGH, IOUtil, and IOHDF5 extensions */
  pughGH = PUGH_pGH (it_info->GH);
  ioUtilGH = (ioGH *) it_info->GH->extensions [CCTK_GHExtensionHandle ("IO")];
  myGH = (ioHDF5GH*) it_info->GH->extensions [CCTK_GHExtensionHandle("IOHDF5")];

  CACTUS_IOHDF5_ERROR (H5Gget_objinfo (group, datasetname, 0, &object_info));
  is_group = object_info.type == H5G_GROUP;
  if (is_group)
    CACTUS_IOHDF5_ERROR (dataset = H5Gopen (group, datasetname));
  else
    CACTUS_IOHDF5_ERROR (dataset = H5Dopen (group, datasetname));

  /* read in the dataset's attributes and verify them */
  if (GetCommonAttributes (it_info->GH, dataset, datasetname,
      it_info->unchunked, &index, &gtype, &timelevel, is_group)) { 
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Ignoring dataset '%s'", datasetname);
    return (0);
  }

  /* if we read in initial data via the file reader interface
     check whether the user wants to have this variable restored */
  if (ioUtilGH->do_inVars && ! ioUtilGH->do_inVars [index]) {
    if (verbose) {
      char *varname = CCTK_FullName (index);

      CCTK_VInfo (CCTK_THORNSTRING, "Ignoring variable '%s' for file reader "
                  "recovery", varname);
      free (varname);
    }
    return (0);
  }

  if (verbose) {
    char *varname = CCTK_FullName (index);

    CCTK_VInfo (CCTK_THORNSTRING, "Restoring variable '%s' (timelevel %d)",
                varname, timelevel);
    free (varname);
  }

  rec_info.it_info = it_info;

  switch (CCTK_VarTypeI (index)) {
    case CCTK_VARIABLE_CHAR:
      rec_info.iohdf5_type = IOHDF5_CHAR;
      rec_info.element_size = sizeof (CCTK_CHAR);
#ifdef CCTK_MPI
      rec_info.mpi_type = PUGH_MPI_CHAR;
#endif
      break;
    case CCTK_VARIABLE_INT:
      rec_info.iohdf5_type = IOHDF5_INT;
      rec_info.element_size = sizeof (CCTK_INT);
#ifdef CCTK_MPI
      rec_info.mpi_type = PUGH_MPI_INT;
#endif
      break;
    case CCTK_VARIABLE_REAL:
      rec_info.iohdf5_type = IOHDF5_REAL;
      rec_info.element_size = sizeof (CCTK_REAL);
#ifdef CCTK_MPI
      rec_info.mpi_type = PUGH_MPI_REAL;
#endif
      break;
    case CCTK_VARIABLE_COMPLEX:
      rec_info.iohdf5_type = myGH->IOHDF5_COMPLEX;
      rec_info.element_size = sizeof (CCTK_COMPLEX);
#ifdef CCTK_MPI
      rec_info.mpi_type = pughGH->PUGH_mpi_complex;
#endif
      break;
    default:
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "Unsupported datatype %d", CCTK_VarTypeI (index));
      return (0);
  }

  /* Read in the data */
  switch (gtype) {
    case CCTK_SCALAR:
      IOHDF5_RestoreGS (dataset, index, timelevel, &rec_info);
      break;
    case CCTK_GF:
    case CCTK_ARRAY:
      IOHDF5_RestoreGA (dataset, index, timelevel, &rec_info);
      break;
    default:
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "Unsupported group type %d", gtype);
      return (1);
  }

  if (is_group)
    CACTUS_IOHDF5_ERROR (H5Gclose (dataset));
  else
    CACTUS_IOHDF5_ERROR (H5Dclose (dataset));

  return (0);
}