CH_HDF5.H

Go to the documentation of this file.

#ifdef CH_LANG_CC
/*
 *      _______              __
 *     / ___/ /  ___  __ _  / /  ___
 *    / /__/ _ \/ _ \/  V \/ _ \/ _ \
 *    \___/_//_/\___/_/_/_/_.__/\___/
 *    Please refer to Copyright.txt, in Chombo's root directory.
 */
#endif

#ifndef _CH_HDF5_H_
#define _CH_HDF5_H_

#define CHOFFSET(object, member) (int)((char*)&(object.member) - (char*)&object)

#ifdef CH_USE_HDF5  // if you don't have CH_USE_HDF5, then this file is useless

#include <iostream>
using std::cout;
using std::endl;

#ifdef CH_MPI
#include <mpi.h>
#endif

#include "LevelData.H"
#include "HDF5Portable.H"
// hdf5 #defines inline.... duh!
#undef inline
#include <string>
#include <map>
#include "RealVect.H"
#include "CH_Timer.H"
#include "LoadBalance.H"
#include "LayoutIterator.H"
#include "Vector.H"
#include "memtrack.H"
#include "FluxBox.H"

#ifdef CH_MULTIDIM
#include "BoxTools_ExternC_Mangler.H" // Generated by lib/utils/multidim/mangle_externs.sh
#endif
#include "NamespaceHeader.H"


template <class T>
void read(T& item, Vector<Vector<char> >& a_allocatedBuffers, const Box& box,
          const Interval& comps)
{
  Vector<void*> b(a_allocatedBuffers.size());
  for (int i=0; i<b.size(); ++i) b[i]= &(a_allocatedBuffers[i][0]);
  read(item, b, box, comps);
}

using std::map;

class HDF5Handle;

// CH_HDF5.H
// ============

/// user-friendly function to write out data on a AMR level
/**
   all data is IN data.
   a_handle:  handle open and ready.
   a_data:    data, what else do you want to know ?
   a_dx:      the grid spacing at this level
   a_dt:      the timestep size that was last completed
   a_time:    the time of this level (might not be the same as other levels)
   a_domain:  the problem domain, represented at this level of refinement
   a_refRatio:the refinement of a_level+1 wrt a_level. for vis systems it
   would probably help if you use 1 for a_level==max_level.
*/
template <class T>
int writeLevel(HDF5Handle& a_handle,
               const int& a_level,
               const T& a_data,
               const Real& a_dx,
               const Real& a_dt,
               const Real& a_time,
               const Box& a_domain,
               const int& a_refRatio,
               const IntVect& outputGhost = IntVect::Zero,
               const Interval& comps = Interval());

template <class T>
int readLevel(HDF5Handle& a_handle,
              const int& a_level,
              LevelData<T>& a_data,
              Real& a_dx,
              Real& a_dt,
              Real& a_time,
              Box& a_domain,
              int& a_refRatio,
              const Interval& a_comps = Interval(),
              bool setGhost = false);

template <class T>
int readLevel(HDF5Handle& a_handle,
              const int& a_level,
              LevelData<T>& a_data,
              RealVect& a_dx,
              Real& a_dt,
              Real& a_time,
              Box& a_domain,
              IntVect& a_refRatio,
              const Interval& a_comps = Interval(),
              bool setGhost = false);

// More basic HDF5 functions.  These can be used at a persons own
// discretion.  Refer to the User's Guide for a description of how these
// functions interact with the convenience functions writeLevel, etc.

/// writes BoxLayout to HDF5 file.
/**
    Writes BoxLayout to HDF5 file.  Only one BoxLayout per group is permitted, this operation overwrites
    previous entries.
    This operation assumes boxes are cell-centered.\\
    returns: success:    0\\
    HDF5 error: negative error code.\\
*/
int write(HDF5Handle& a_handle,
          const BoxLayout& a_layout,
          const std::string& name = "boxes");

/// writes a BoxLayoutData<T> to an HDF5 file.
/**
   writes a BoxLayoutData<T> to an HDF5 file.\\
   returns: success:    0\\
   HDF5 error: negative error code.\\
*/
template <class T>
int write(HDF5Handle& a_handle,
          const BoxLayoutData<T>& a_data,
          const std::string& a_name,
          IntVect outputGhost = IntVect::Zero,
          const Interval& comps = Interval(),
          bool  newForm = false);

/// writes a LevelData<T> to an HDF5 file.
/**
   Writes a LevelData<T> to an HDF5 file.
   the DisjointBoxLayout is not written out with the data, the user is required to
   handle that object seperately. (see the "read" interface below).\\
   returns: success:    0\\
   HDF5 error: negative error code.\\
*/
template <class T>
int write(HDF5Handle& a_handle,
          const LevelData<T>& a_data,
          const std::string& a_name,
          const IntVect& outputGhost = IntVect::Zero,
          const Interval& comps = Interval());

/// reads Vector<Box> from location specified by a_handle.
/**
     Reads BoxLayout from the group specified by a_handle.
     Only one BoxLayout per group is permitted, this operation overwrites.
     This operation assumes boxes are cell-centered.\\
     returns: success:    0\\
     HDF5 error: negative error code.\\

     Arg name refers to the name of an HDF5 dataset in which the boxes you want
     are stored in a particular format, which is the one used to output DataLayouts
     to HDF5.  Here is an example of that format (as shown by h5dump):

      DATASET "datalayout" {
         DATATYPE  H5T_COMPOUND {
            H5T_STD_I32LE "lo_i";
            H5T_STD_I32LE "lo_j";
            H5T_STD_I32LE "hi_i";
            H5T_STD_I32LE "hi_j";
         }
         DATASPACE  SIMPLE { ( 2 ) / ( 2 ) }
         DATA {
         (0): {
               0,
               0,
               1,
               0
            },
         (1): {
               2,
               0,
               2,
               1
            }
         }
      }

*/
int read(HDF5Handle& a_handle,
         Vector<Box>& boxes,
         const std::string& name = "boxes");

/// reads the set of Boxes out from the level_* groups of a Chombo HDF5 AMR file
/**
   goes to all groups named level_n for level n = 0 to numLevel and fills the
   vector of boxes.  CH_HDF5.Handle must be set to the root group (the default location when the
   file is just opened.

   A handy *get everything* version of int read(HDF5Handle& a_handle, Vector<Box>& boxes)
*/
int readBoxes(HDF5Handle& a_handle,
              Vector<Vector<Box> >& boxes);

/// FArrayBox-at-a-time read function.  FArrayBox gets redefined in the function.  Reads data field named by a_dataName.
/**  FArrayBox gets redefined in the function.
     it will be sized to the box at the [level,boxNumber] position and have components
     the go from [0,a_components.size()-1]

     some meta-data has to be read again and again internally to do this function, but
     it should make some out-of-core tools possible.
*/

int readFArrayBox(HDF5Handle& a_handle,
                  FArrayBox&  a_fab,
                  int a_level,
                  int a_boxNumber,
                  const Interval& a_components,
                  const std::string& a_dataName = "data" );

/// read BoxLayoutData named a_name from location specified by a_handle.
/**
    Read BoxLayoutData named a_name from location specified by a_handle.  User must supply the correct BoxLayout for this function if redefineData == true.  \\
    returns: success:      0\\
    bad location: 1\\
    HDF5 error:   negative error code.\\
*/
template <class T>
int read(HDF5Handle& a_handle,
         BoxLayoutData<T>& a_data,
         const std::string& a_name,
         const BoxLayout& a_layout,
         const Interval&  a_comps = Interval(),
         bool redefineData = true);

/// read LevelData named a_name from location specified by a_handle.
/**

    Read LevelData named a_name from location specified by a_handle.
    User must supply the correct BoxLayout for this function.

    Arg a_name is significant: the HDF5 group to which a_handle is set is
    is assumed to contain a dataset called <a_name>:datatype=<some integer>,
    a dataset called <a_name>:offsets=<some integer>, and a subgroup named
    <a_name>_attributes.  You will have all these items if you dumped your
    LevelData out using the corresponding write() function defined here.

    If arg redefineData==false, then the user must pass in a valid  LevelData.
    Otherwise, this function figures out how many components and ghost cells there
    are, and allocates the correct amount of space.  The actual FArray data held
    by the LevelData gets filled in here, regardless of redefineData; "redefine"
    alludes to the family of define() functions.

    returns: success:      0\\
    bad location: 1\\
    HDF5 error:   negative error code.\\
*/
template <class T>
int read(HDF5Handle& a_handle,
         LevelData<T>& a_data,
         const std::string& a_name,
         const DisjointBoxLayout& a_layout,
         const Interval& a_comps = Interval(),
         bool redefineData = true);

/// Handle to a particular group in an HDF file.
/**
    HDF5Handle is a handle to a particular group in an HDF file.  Upon
    construction, group is defined to be the root.  All data is
    written and read assuming the native representations for data on
    the architecture it is running on.  When a file is opened, these
    settings are checked and an error is flagged when things don't
    match up.  It is the USER'S responsibility to close() this object
    when it is no longer needed.

*/
class HDF5Handle
{
public:
  ///
  /**
     Enumeration of opening modes for HDF5 files.  \\

     CREATE: file is created if it didn't exist, or an existing file of
     the same name is clobbered.\\

     CREATE_SERIAL: in serial execution this is equivalent to CREATE. In parallel this opens
     a file just on this particular calling processor.  Working with LevelData and this kind of
     file is difficult to get right and most users will not use this mode.  \\

     OPEN_RDONLY: existing file is opened in read-only mode.  If the
     file doesn't already exist then open fails and isOpen() returns
     false.\\

     OPEN_RDWR: existing file is opened in read-write mode.  If the file
     doesn't already exist then open fails and isOpen() returns false.\\

  */
  enum mode
  {
    CREATE,
    CREATE_SERIAL,
    OPEN_RDONLY,
    OPEN_RDWR
  };

  ///    {\bf constructor}

  ///
  /**
     Default constructor.  User must call open() prior to using
     constructed object.
  */
  HDF5Handle();

  ///
  /** Opens file and sets the current group to the root "/" group. \\

  if mode == CREATE, then file is created if it didn't exist, or an
  existing file of the same name is clobbered.\\

  if mode == OPEN_*, then existing file is opened, if the file doesn't
  already exist then open fails and isOpen() returns false.\\

  Writes basic information such as SpaceDim and testReal to a global
  group. The global group is named Chombo_global in the call signature
  where no globalGroupName is passed.\\
  */
  HDF5Handle(
        const std::string& a_filename,
        mode a_mode,
        const char *a_globalGroupName="Chombo_global");

  // HDF5Handle(const std::string& a_filename, mode a_mode);

  ~HDF5Handle();

  /// {\bf File functions}

  ///
  /**
      Opens file and sets the current group of this HDF5Handle to the
      root "/" group.  File that this HDF5Handle previously pointed at is
      NOT closed, that is the users responsibility.\\

      if mode == CREATE, then file is created if it didn't exist, or an
      existing file of the same name is clobbered.\\

      if mode == OPEN_*, then existing file is opened, if the file doesn't
      already exist then open fails and isOpen() returns false.\\

      Writes basic information such as SpaceDim and testReal to a global
      group. The global group is named Chombo_global in the call signature
      where no globalGroupName is passed.\\

      returns:\\
      0  on success\\
      negative number if file open failed (return code from HDF5)\\
      1  file does not appear to contain datacheck info, probably not a Chombo file\\
      2  on data bit size differences between code and file.\\

      aborts on SpaceDim not matching between code and file\\

  */
  int open(
        const std::string& a_filename,
        mode a_mode,
        const char *a_globalGroupName="Chombo_global");

  // int open(const std::string& a_filename, mode a_mode);

  ///
  /**
     A NULL or failed constructed HDF5Handle will return false.
  */
  bool isOpen() const;

  ///
  /**
     Closes the file.  Must be called to close file.  Files are not
     automatically closed.

  */
  void close();

  /// {\bf Group functions}

  ///
  /**
     Sets the current group to be "/level_x" where x=a_level.
  */
  void setGroupToLevel(int a_level);

  ///
  /**
      Set group to users choice, referenced from file root.
      groupAbsPath will look like a Unix file path:
      "/mySpecialData/group1/" "/" is the root group of the
      file. returns a negative value on failure

  */
  int setGroup(const std::string& groupAbsPath);

  ///
  /**
      Add the indicated string to the group path.  For example, if
      getGroup() returns "/foo" then after pushGroup("bar"), getGroup()
      will return "/foo/bar".
      Return value is whatever the internal setGroup returns.
  */
  int pushGroup( const std::string& grp );

  ///
  /**
      Pop off the last element of the group path.  For example, "/foo/bar"
      becomes "/foo".  It's an error to call this function if the group,
      going in, is "/".
      Return value is whatever the internal setGroup returns, when we call
      it to reset the group's path.
  */
  int popGroup();

  ///
  /**
     Returns name of current group.  groupAbsPath will look like a Unix
     file path: "/mySpecialData/group1/" "/" is the root group of the
     file.

  */
  const std::string& getGroup() const;

  const hid_t& fileID() const;
  const hid_t& groupID() const;
  static hid_t box_id;
  static hid_t intvect_id;
  static hid_t realvect_id;
  static map<std::string, std::string> groups;

private:

  HDF5Handle(const HDF5Handle&);
  HDF5Handle& operator=(const HDF5Handle&);

  hid_t         m_fileID;
  hid_t         m_currentGroupID;
  bool          m_isOpen;
  std::string   m_filename; // keep around for debugging
  std::string   m_group;
  int           m_level;

  //  static hid_t  file_access;
  static bool   initialized;
  static void   initialize();

};

/// data to be added to HDF5 files.
/**
   HDF5HeaderData is a wrapper for some data maps to be added to HDF5
   files.  instead of an overdose of access functions, the maps are
   made public and they can be manipulated by the user at will.  They
   maintain type safety.

   to add a Real data entry, a user can simply program as follows:

   <PRE>
   Real dx;
   .
   .
   HDF5HeaderData metaData;
   metaData.m_real["dx"] = dx;
   </PRE>

   If "dx" already existed, then it is overwritten, otherwise an entry is
   created and added with the new value;

   To search for entries, the user does the following:
   <PRE>

   HDF5HeaderData metaData;
   HDF5Handle currentStep(filename);
   currentStep.setGroupToLevel(0);
   metaData.readFromFile(currentStep);
   if (metaData.m_intvect.find("ghost") != metaData.m_intvect.end())
   ghost = metaData.m_intvect["ghost"];
   else
   ghost = defaultGhostIntVect;

   </PRE>

   A user can skip the check for existence if they have reason to "know" the
   data will be there.  It is just good coding practice.

   To erase an entry, you can use:
   <PRE>
   metaData.m_real.erase("dx");
   </PRE>

*/
class HDF5HeaderData
{
public:

  ///
  /**
      Writes this HDF5HeaderData's current attribute list to the
      current group in 'file.'  Returns 0 on success, returns the
      error code from HDF5 on failure.

  */
  int writeToFile(HDF5Handle& file) const;

  ///
  /**
      Reads into this HDF5HeaderData's attribute list from file.  Read
      process is add/change, does not remove key-value pairs. Reads
      from current group.  Returns 0 on success, positive number if a
      particular member of group caused an error, negative on general
      error.

  */
  int readFromFile(HDF5Handle& file);

  ///
  void clear();

  ///
  map<std::string, Real>        m_real;

  ///
  map<std::string, int>         m_int;

  ///
  map<std::string, std::string> m_string;

  ///
  map<std::string, IntVect>     m_intvect;

  ///
  map<std::string, Box>         m_box;

  ///
  map<std::string, RealVect>    m_realvect;

  //users should not need these functions in general

  int writeToLocation(hid_t loc_id) const;
  int readFromLocation(hid_t loc_id);

  /// useful for debugging.  dumps contents to std::cout
  void dump() const;

private:
  static herr_t attributeScan(hid_t loc_id, const char *name, void *opdata);
};

extern "C"
{
  herr_t HDF5HeaderDataattributeScan(hid_t loc_id, const char *name, void *opdata);
}

std::ostream& operator<<(std::ostream& os, const HDF5HeaderData& data);

//=============================================================================
//
// end of declarations.
//
//=============================================================================

#if ( H5_VERS_MAJOR == 1 && H5_VERS_MINOR > 6 )
typedef hsize_t ch_offset_t;
#else
#if ( H5_VERS_MAJOR == 1 && H5_VERS_MINOR == 6 && H5_VERS_RELEASE >= 4 )
typedef hsize_t ch_offset_t;
#else
typedef hssize_t ch_offset_t;
#endif
#endif


template<class T>
hid_t H5Type(const T* dummy);

template< >
hid_t H5Type(const int* dummy);

template< >
hid_t H5Type(const long long* dummy);

template< >
hid_t H5Type(const float* dummy);

template< >
hid_t H5Type(const double* dummy);

template< >
hid_t H5Type(const Box* dummy);

template< >
hid_t H5Type(const RealVect* dummy);

template< >
hid_t H5Type(const IntVect* dummy);

template<class T>
hid_t H5Type(const T* dummy)
{
  // no such definition;
  MayDay::Error(" H5Type(const T* dummy)");
  return -4;
}


void createData(hid_t& a_dataset,
                hid_t& a_dataspace,
                HDF5Handle& handle,
                const std::string& name,
                hid_t type,
                hsize_t size);

template <class T>
void createDataset(hid_t& a_dataset,
                   hid_t& a_dataspace,
                   HDF5Handle& handle,
                   const std::string& name,
                   const T* dummy,
                   hsize_t size)
{
  createData(a_dataset, a_dataspace, handle, name, H5Type(dummy), size);

}

void writeDataset(hid_t a_dataset,
                  hid_t a_dataspace,
                  const void* start,
                  ch_offset_t off,
                  hsize_t  count);

void readDataset(hid_t a_dataset,
                 hid_t a_dataspace,
                 void* start,
                 ch_offset_t off,
                 hsize_t  count);

// non-user code used in implementation of communication

struct OffsetBuffer
{
  Vector<int> index;
  Vector<Vector<int> > offsets;
  void operator=(const OffsetBuffer& rhs);
};

ostream& operator<<(ostream& os, const OffsetBuffer& ob);


#include "NamespaceFooter.H"

#include "BaseNamespaceHeader.H"

#include "NamespaceVar.H"

//OffsetBuffer specialization of linearSize
template < >
int linearSize(const CH_XDIR::OffsetBuffer& a_input);

//OffsetBuffer specialization of linearIn
template < >
void linearIn(CH_XDIR::OffsetBuffer& a_outputT, const void* const a_inBuf);

//OffsetBuffer specialization of linearOut
template < >
void linearOut(void* const a_outBuf, const CH_XDIR::OffsetBuffer& a_inputT);

template < > int linearSize(const Vector<CH_XDIR::OffsetBuffer>& a_input);
template < > void linearIn(Vector<CH_XDIR::OffsetBuffer>& a_outputT, const void* const inBuf);
template < > void linearOut(void* const a_outBuf, const Vector<CH_XDIR::OffsetBuffer>& a_inputT);


#include "BaseNamespaceFooter.H"
#include "NamespaceHeader.H"

// First, template specializations for read/write for FArrayBox.

template <>
inline void dataTypes(Vector<hid_t>& a_types, const BaseFab<int>& dummy)
{
  a_types.resize(1);
  a_types[0] = H5T_NATIVE_INT;
}

template <>
inline void dataTypes(Vector<hid_t>& a_types, const BaseFab<char>& dummy)
{
  a_types.resize(1);
  a_types[0] = H5T_NATIVE_CHAR;
}

/* since many compilers choke on the proper template specialization syntax
   I am forced to pass in a dummy specialization argument
*/
template <>
inline void dataTypes(Vector<hid_t>& a_types, const FArrayBox& dummy)
{
  a_types.resize(1);
  a_types[0] = H5T_NATIVE_REAL;
}

template <>
inline void dataTypes(Vector<hid_t>& a_types, const FluxBox& dummy)
{
  a_types.resize(1);
  a_types[0] = H5T_NATIVE_REAL;
}

template <>
inline void dataSize(const BaseFab<int>& item, Vector<int>& a_sizes,
                     const Box& box, const Interval& comps)
{
  a_sizes[0] = box.numPts() * comps.size();
}

template <>
inline void dataSize(const BaseFab<char>& item, Vector<int>& a_sizes,
                     const Box& box, const Interval& comps)
{
  a_sizes[0] = box.numPts() * comps.size();
}

template <>
inline void dataSize(const FArrayBox& item, Vector<int>& a_sizes,
                     const Box& box, const Interval& comps)
{
  a_sizes[0] = box.numPts() * comps.size();
}

template <>
inline void dataSize(const FluxBox& item, Vector<int>& a_sizes,
                     const Box& box, const Interval& comps)
{
  int& t = a_sizes[0];
  t = 0;
  for (int dir =0; dir<CH_SPACEDIM; dir++)
    {
      Box edgeBox(surroundingNodes(box,dir));
      t += edgeBox.numPts()*comps.size();
    }
}

template <>
inline const char* name(const FArrayBox& a_dummySpecializationArg)
{
  // Attempt to get rid of warnings on IBM...
  //static const char* name = "FArrayBox";
  const char* name = "FArrayBox";
  return name;
}

template <>
inline const char* name(const BaseFab<int>& a_dummySpecializationArg)
{
  // Attempt to get rid of warnings on IBM...
  //static const char* name = "BaseFab<int>";
  const char* name = "BaseFab<int>";
  return name;
}

template <>
inline const char* name(const BaseFab<char>& a_dummySpecializationArg)
{
  // Attempt to get rid of warnings on IBM...
  //static const char* name = "BaseFab<int>";
  const char* name = "BaseFab<char>";
  return name;
}
//
// now, generic, non-binary portable version of template functions
// for people who just want ot rely on linearIn/linearOut

template <class T>
inline void dataTypes(Vector<hid_t>& a_types, const T& dummy)
{
  a_types.resize(1);
  a_types[0] = H5T_NATIVE_CHAR;
}

template <class T>
inline void dataSize(const T& item, Vector<int>& a_sizes,
                     const Box& box, const Interval& comps)
{
  a_sizes[0] = item.size(box, comps);
}

template <class T>
inline void write(const T& item, Vector<void*>& a_allocatedBuffers,
                  const Box& box, const Interval& comps)
{
  item.linearOut(a_allocatedBuffers[0], box, comps);
}

template <class T>
inline void read(T& item, Vector<void*>& a_allocatedBuffers,
                 const Box& box, const Interval& comps)
{
  item.linearIn(a_allocatedBuffers[0], box, comps);
}

template <class T>
inline const char* name(const T& a_dummySpecializationArg)
{
  static const char* name = "unknown";
  return name;
}

template <class T>
void getOffsets(Vector<Vector<long long> >& offsets, const BoxLayoutData<T>& a_data,
                int types, const Interval& comps, const IntVect& outputGhost)
{
  CH_TIME("getOffsets");
  const BoxLayout& layout = a_data.boxLayout();
  offsets.resize(types, Vector<long long>(layout.size()+1));
  //  offsets.resize(layout.size() + 1, Vector<long long>(types));
  for (int t=0; t<types; t++) offsets[t][0] = 0;
  Vector<int> thisSize(types);
  if (T::preAllocatable()==0)
    { // static preAllocatable
      T dummy;
      unsigned int index = 1;
      for (LayoutIterator it(layout.layoutIterator()); it.ok(); ++it)
        {
          Box region = layout[it()];
          region.grow(outputGhost);
          dataSize(dummy, thisSize, region, comps);
          for (unsigned int i=0; i<thisSize.size(); ++i)
            {
              //offsets[index][i] = offsets[index-1][i] + thisSize[i];
              offsets[i][index] = offsets[i][index-1] + thisSize[i];
            }
          ++index;
        }
    }
  else
    { // symmetric and dynamic preallocatable need two pass I/O
      OffsetBuffer buff;
      //int index = 0;
      for (DataIterator dit(a_data.dataIterator()); dit.ok(); ++dit)
        {
          int index = a_data.boxLayout().index(dit());
          //int index = dit().intCode();
          buff.index.push_back(index);
          Box region = layout[dit()];
          region.grow(outputGhost);
          dataSize(a_data[dit()], thisSize, region, comps);
          buff.offsets.push_back(thisSize);
        }
      Vector<OffsetBuffer> gathering(numProc());
      gather(gathering, buff, uniqueProc(SerialTask::compute));
      broadcast(gathering,  uniqueProc(SerialTask::compute));
      //  pout() << gathering<<endl;
      for (int i=0; i<numProc(); ++i)
        {
          OffsetBuffer& offbuf = gathering[i];
          for (int num=0; num<offbuf.index.size(); num++)
            {
              int index = offbuf.index[num];
              for (unsigned int j=0; j<types; ++j)
                {
                  //offsets[index+1][j] = offbuf.offsets[num][j];
                  offsets[j][index+1] = offbuf.offsets[num][j];
                }
            }
        }
      for (int i=0; i<layout.size(); i++)
        {
          for (unsigned int j=0; j<types; ++j)
            {
              //offsets[i+1][j] += offsets[i][j];
              offsets[j][i+1] += offsets[j][i];
            }
        }
    }

  // pout() << offsets<<endl;
}

//==================================================================
//
// Now, linear IO routines for a BoxLayoutData of T
//
template <class T>
int write(HDF5Handle& a_handle, const BoxLayoutData<T>& a_data,
          const std::string& a_name, IntVect outputGhost,
          const Interval& in_comps, bool newForm)
{
  CH_TIME("write_Level");
  int ret = 0;

  Interval comps(in_comps);
  if ( comps.size() == 0) comps = a_data.interval();
  T dummy; // used for preallocatable methods for dumb compilers.
  Vector<hid_t> types;
  dataTypes(types, dummy);

  Vector<Vector<long long> > offsets;
  Vector<long long> bufferCapacity(types.size(), 1);  // noel (was 0)
  Vector<void*> buffers(types.size(), NULL);

  getOffsets(offsets, a_data, types.size(), comps, outputGhost);

  // create datasets collectively.
  hsize_t flatdims[1];
  char dataname[100];
  Vector<hid_t> dataspace(types.size());
  Vector<hid_t> dataset(types.size());

  herr_t err;
  hsize_t count[1];
  ch_offset_t offset[1];
  CH_assert(!(newForm && types.size() != 1));

  for (unsigned int i=0; i<types.size(); ++i)
    {
      flatdims[0] = offsets[i][offsets[i].size()-1];
      if (newForm)
      {
        if (a_name == "M")
        {
          strcpy(dataname, "Mask");
        }
        else
        {
          sprintf(dataname, "%sRegular", a_name.c_str());
        }
      }
      else
      {
        sprintf(dataname, "%s:datatype=%i",a_name.c_str(), i);
      }
      {
        CH_TIME("H5Screate");
        dataspace[i]      = H5Screate_simple(1, flatdims, NULL);
      }
      CH_assert(dataspace[i] >=0);
      {
        CH_TIME("H5Dcreate");
        dataset[i]        = H5Dcreate(a_handle.groupID(), dataname,
                                      types[i],
                                      dataspace[i], H5P_DEFAULT);
      }
      CH_assert(dataset[i] >= 0);
    }

  hid_t offsetspace, offsetData;
  for (unsigned int i=0; i<types.size(); ++i)
    {
      flatdims[0] =  offsets[i].size();
      if (newForm)
      {
        if (a_name == "M")
        {
          strcpy(dataname, "MaskOffsets");
        }
        else
        {
          sprintf(dataname, "%sOffsets",a_name.c_str());
        }
      }
      else
      {
        sprintf(dataname, "%s:offsets=%i",a_name.c_str(), i);
      }
      {
        CH_TIME("H5S_H5D_offsets_create");
        offsetspace =  H5Screate_simple(1, flatdims, NULL);
        CH_assert(offsetspace >= 0);
        offsetData  =   H5Dcreate(a_handle.groupID(), dataname,
                                  H5T_NATIVE_LLONG, offsetspace, H5P_DEFAULT);
        CH_assert(offsetData >= 0);
      }
      if (procID() == 0)
        {
          CH_TIME("WriteOffsets");
          hid_t memdataspace = H5Screate_simple(1, flatdims, NULL);
          CH_assert(memdataspace >= 0);
          err = H5Dwrite(offsetData, H5T_NATIVE_LLONG, memdataspace, offsetspace,
                         H5P_DEFAULT, &(offsets[i][0]));
          CH_assert(err >= 0);
          H5Sclose(memdataspace);
        }
      {
        CH_TIME("H5S_H5D_offsets_close");
        H5Sclose(offsetspace);
        H5Dclose(offsetData);
      }
    }

  // write BoxLayoutData attributes into Dataset[0]
  if (!newForm)
  {
    CH_TIME("WriteAttributes");
    HDF5HeaderData info;
    info.m_int["comps"] = comps.size();
    info.m_string["objectType"] = name(dummy);
    std::string group = a_handle.getGroup();
    a_handle.setGroup(group+"/"+a_name+"_attributes");
    info.writeToFile(a_handle);
    a_handle.setGroup(group);
  }

  // collective operations finished, now perform parallel writes
  // to specified hyperslabs.

  Vector<size_t> type_size(types.size());
  for (unsigned int i=0; i<types.size(); ++i)
    {
      type_size[i] = H5Tget_size(types[i]);
    }

  Vector<int> thisSize(types.size());

  // step 1, create buffer big enough to hold the biggest linearized T
  // that I will have to output.

  //  pout()<<"offsets ";
  for (int i=0; i<offsets[0].size(); i++)
  {
    //    pout()<<" "<<offsets[0][i];
  }
  // pout()<<"\n";
  {
    CH_TIME("ComputeBufferCapacity");
    for (DataIterator it = a_data.dataIterator(); it.ok(); ++it)
      {
        unsigned int index = a_data.boxLayout().index(it());
        for (unsigned int i=0; i<types.size(); ++i)
          {
            long long size = (offsets[i][index+1] - offsets[i][index])
              * type_size[i];
            // pout()<<"index:size "<<index<<" "<<size<<"\n";
            CH_assert(size >= 0);
            if (size > bufferCapacity[i]) // grow buffer if necessary.....
              {
                bufferCapacity[i] = size;
              }
          }
      }
  }
  // CH_assert(bufferCapacity[0] > 1);
  for (unsigned int i=0; i<types.size(); ++i)
    {
      CH_TIME("mallocMT_buffers");
      buffers[i] = mallocMT(bufferCapacity[i]);
      if (buffers[i] == NULL)
      {
        pout() << " i=" << i
               << " types.size() = " << types.size()
               << " bufferCapacity[i] = " << (int)bufferCapacity[i]
               << endl;
        MayDay::Error("memory error in buffer allocation write");
      }
    }

  // Step 2.  actually a) write each of my T objects into the
  // buffer, then b) write that buffered data out to the
  // write position in the data file using hdf5 hyperslab functions.
#ifdef TRY_MPI_COLLECTIVES_
  hid_t DXPL = H5Pcreate(H5P_DATASET_XFER);
  H5Pset_dxpl_mpio(DXPL, H5FD_MPIO_COLLECTIVE);
#endif
  {
    CH_TIME("linearize_H5Dwrite");
    for (DataIterator it = a_data.dataIterator(); it.ok(); ++it)
      {
        const T& data = a_data[it()];
        unsigned int index = a_data.boxLayout().index(it());
        Box box = a_data.box(it());
        box.grow(outputGhost);
        {       
          CH_TIMELEAF("linearize");
          write(data, buffers, box, comps); //write T to buffer
        }
        for (unsigned int i=0; i<types.size(); ++i)
          {
            offset[0] = offsets[i][index];
            count[0] = offsets[i][index+1] - offset[0];
            hid_t memdataspace;
            if (count[0] > 0)
              {
              err =  H5Sselect_hyperslab(dataspace[i], H5S_SELECT_SET,
                                         offset, NULL,
                                         count, NULL);
              CH_assert(err >= 0);
              memdataspace = H5Screate_simple(1, count, NULL);
              CH_assert(memdataspace >= 0);
            }
          else
            {
              H5Sselect_none(dataspace[i]);
              H5Sselect_none(memdataspace);
            }
          {
            CH_TIMELEAF("H5Dwrite");
#ifdef TRY_MPI_COLLECTIVES_
            err = H5Dwrite(dataset[i], types[i], memdataspace, dataspace[i],
                           DXPL, buffers[i]);
#else
            err = H5Dwrite(dataset[i], types[i], memdataspace, dataspace[i],
                           H5P_DEFAULT, buffers[i]);
#endif
          }
          CH_assert(err >= 0);
          H5Sclose(memdataspace);
          if (err < 0)
            {
              ret = err;
              goto cleanup;
            }
          
        }
    }
  }
#ifdef TRY_MPI_COLLECTIVES_
  H5Pclose(DXPL);
#endif
  
  // OK, clean up data structures
  
 cleanup:
  for (unsigned int i=0; i<types.size(); ++i)
    {
      {
        CH_TIME("freeMT");
        freeMT(buffers[i]);
      }
      {
        CH_TIME("H5Sclose");
        H5Sclose(dataspace[i]);
      }
      {
        CH_TIME("H5Dclose");
        H5Dclose(dataset[i]);
      }
    }
  return ret;
  
}

template <class T>
int write(HDF5Handle& a_handle, const LevelData<T>& a_data,
          const std::string& a_name, const IntVect& outputGhost, const Interval& in_comps)
{
  HDF5HeaderData info;
  info.m_intvect["ghost"] = a_data.ghostVect();
  IntVect og(outputGhost);
  og.min(a_data.ghostVect());
  info.m_intvect["outputGhost"] = og;
  std::string group = a_handle.getGroup();
  a_handle.setGroup(group+"/"+a_name+"_attributes");
  info.writeToFile(a_handle);
  a_handle.setGroup(group);
  return write(a_handle, (const BoxLayoutData<T>&)a_data, a_name, og, in_comps);
}

template <class T>
int read(HDF5Handle& a_handle, LevelData<T>& a_data, const std::string& a_name,
         const DisjointBoxLayout& a_layout, const Interval& a_comps, bool a_redefineData)
{
  if (a_redefineData)
    {
      HDF5HeaderData info;
      std::string group = a_handle.getGroup();
      if (a_handle.setGroup(group+"/"+a_name+"_attributes"))
        {
          std::string message = "error opening "
                                +a_handle.getGroup()+"/"+a_name+"_attributes" ;
          MayDay::Warning(message.c_str());
          return 1;
        }
      info.readFromFile(a_handle);
      a_handle.setGroup(group);
      int ncomp =  info.m_int["comps"];
      IntVect ghost = info.m_intvect["ghost"];
      if (a_comps.end() > 0 && ncomp < a_comps.end())
        {
          MayDay::Error("attempt to read component interval that is not available");
        }
      if (a_comps.size() == 0)
        a_data.define(a_layout, ncomp, ghost);
      else
        a_data.define(a_layout, a_comps.size(), ghost);
    }
  return read(a_handle, (BoxLayoutData<T>&)a_data, a_name, a_layout, a_comps, false);

}

template <class T>
int read(HDF5Handle& a_handle, BoxLayoutData<T>& a_data, const std::string& a_name,
         const BoxLayout& a_layout, const Interval& a_comps, bool a_redefineData)
{
  CH_TIME("read");
  int ret = 0; // return value;

  herr_t err;

  char dataname[100];
  hsize_t count[1];
  ch_offset_t offset[1];
  Vector<Vector<long long> > offsets;

  T dummy;
  Vector<hid_t> types;
  dataTypes(types, dummy);
  Vector<hid_t> dataspace(types.size());
  Vector<hid_t> dataset(types.size());
  offsets.resize(types.size(), Vector<long long>(a_layout.size() +1));

  //Vector<int> bufferCapacity(types.size(), 500);
  //Vector<void*> buffers(types.size(), NULL);
  Vector<Vector<char> > buffers(types.size(), 500);

  for (unsigned int i=0; i<types.size(); ++i)
    {
      sprintf(dataname, "%s:datatype=%i",a_name.c_str(), i);
      dataset[i]        = H5Dopen(a_handle.groupID(), dataname);
      if (dataset[i] < 0)
      {
        MayDay::Warning("dataset open failure"); return dataset[i];
      }
      dataspace[i]      = H5Dget_space(dataset[i]);
      if (dataspace[i] < 0)
      {
        MayDay::Warning("dataspace open failure"); return dataspace[i];
      }
    }

  hid_t offsetspace, offsetData;
  hsize_t flatdims[1];
  for (unsigned int i=0; i<types.size(); ++i)
    {
      flatdims[0] =  offsets[i].size();
      sprintf(dataname, "%s:offsets=%i",a_name.c_str(), i);
      offsetspace =  H5Screate_simple(1, flatdims, NULL);
      CH_assert(offsetspace >= 0);
      offsetData  =   H5Dopen(a_handle.groupID(), dataname);
      CH_assert(offsetData >= 0);
      hid_t memdataspace = H5Screate_simple(1, flatdims, NULL);
      CH_assert(memdataspace >= 0);
      err = H5Dread(offsetData, H5T_NATIVE_LLONG, memdataspace, offsetspace,
                    H5P_DEFAULT, &(offsets[i][0]));
      CH_assert(err >=0);
      H5Sclose(memdataspace);
      H5Sclose(offsetspace);
      H5Dclose(offsetData);
    }

  HDF5HeaderData info;
  std::string group = a_handle.getGroup();
  if (a_handle.setGroup(a_handle.getGroup()+"/"+a_name+"_attributes"))
    {
      std::string message = "error opening "+a_handle.getGroup()+"/"+a_name ;
      MayDay::Warning(message.c_str());
      return 1;
    }

  info.readFromFile(a_handle);
  a_handle.setGroup(group);
  int ncomps = info.m_int["comps"];
  IntVect outputGhost(IntVect::Zero); // backwards file compatible mode.
  if (info.m_intvect.find("outputGhost") != info.m_intvect.end())
    {
      outputGhost = info.m_intvect["outputGhost"];
    }
  if (ncomps <= 0)
  {
    MayDay::Warning("ncomps <= 0 in read");
    return ncomps;
  }

  if (a_redefineData)
  {
    if (a_comps.size() != 0)
    {
      a_data.define(a_layout, a_comps.size());
    }
    else
    {
      a_data.define(a_layout, ncomps);
    }
  }

  Interval comps(0, ncomps-1);
  //huh?
  //  if (a_comps.size() != 0)  comps = Interval(0, a_comps.size());

  //  getOffsets(offsets, a_data, types.size(), comps);

  Vector<size_t> type_size(types.size());
  for (unsigned int i=0; i<types.size(); ++i)
    {
      type_size[i] = H5Tget_size(types[i]);

    }

  Vector<int> thisSize(types.size());
  DataIterator it = a_data.dataIterator();
  for( ; it.ok(); ++it)
    {
      CH_TIMELEAF("H5Dread");
      T& data = a_data[it()];
      unsigned int index = a_data.boxLayout().index(it());
      Box box = a_data.box(it());

      for (unsigned int i=0; i<types.size(); ++i)
        {

          offset[0] = offsets[i][index];
          count[0] = offsets[i][index+1] - offset[0];
          if (count[0] > 0)
          {
            int size = count[0] * type_size[i];
            while (size > buffers[i].size())
              {
                buffers[i].resize(2*buffers[i].size());
            }
            
            err =  H5Sselect_hyperslab(dataspace[i], H5S_SELECT_SET,
                                       offset, NULL,
                                       count, NULL);
            CH_assert(err >= 0);
            hid_t memdataspace = H5Screate_simple(1, count, NULL);
            CH_assert(memdataspace >= 0);
            err = H5Dread(dataset[i], types[i], memdataspace, dataspace[i],
                          H5P_DEFAULT, &(buffers[i][0]));
            CH_assert(err >= 0);
            H5Sclose(memdataspace);
            if (err < 0)
              {
              ret = err;
              goto cleanup;
              }
          }
        }
      box.grow(outputGhost);
      read(data, buffers,  box, comps);
    }
//   if(it.size()==0)
//     {
//       //  try doing an empty H5Dread operation to make the collective system happier.
//       for (unsigned int i=0; i<types.size(); ++i)
//         {
//        hid_t filespace, memspace;
//        H5Sselect_none(filespace);
//        H5Sselect_none(memspace);
//        err = H5Dread(dataset[i], types[i], memspace, filespace,
//                      H5P_DEFAULT, 0);
//        H5Sclose(filespace);
//        H5Sclose(memspace);
//      }
//     }
                        
 cleanup:
  for (unsigned int i=0; i<types.size(); ++i)
    {
      // freeMT(buffers[i]);
      H5Sclose(dataspace[i]);
      H5Dclose(dataset[i]);
    }
  return ret;
}

template <class T>
int writeLevel(HDF5Handle& a_handle,
               const int&  a_level,
               const T& a_data,
               const Real& a_dx,
               const Real& a_dt,
               const Real& a_time,
               const Box&  a_domain,
               const int&  a_refRatio,
               const IntVect& outputGhost,
               const Interval& comps)
{
  int error;
  char levelName[10];
  std::string currentGroup = a_handle.getGroup();
  sprintf(levelName, "/level_%i",a_level);
  error = a_handle.setGroup(currentGroup + levelName);
  if (error != 0) return 1;

  HDF5HeaderData meta;
  meta.m_real["dx"] = a_dx;
  meta.m_real["dt"] = a_dt;
  meta.m_real["time"] = a_time;
  meta.m_box["prob_domain"] = a_domain;
  meta.m_int["ref_ratio"] = a_refRatio;

  error = meta.writeToFile(a_handle);
  if (error != 0) return 2;

  error = write(a_handle, a_data.boxLayout());
  if (error != 0) return 3;

  error = write(a_handle, a_data, "data", outputGhost, comps);
  if (error != 0) return 4;

  a_handle.setGroup(currentGroup);

  return 0;
}


template <class T>
int writeLevel(HDF5Handle& a_handle,
               const int&  a_level,
               const T& a_data,
               const RealVect& a_dx, // dx for each direction
               const Real& a_dt,
               const Real& a_time,
               const Box&  a_domain,
               const IntVect&  a_refRatios, // ref ratio for each direction
               const IntVect& outputGhost,
               const Interval& comps)
{
  int error;
  char levelName[10];
  std::string currentGroup = a_handle.getGroup();
  sprintf(levelName, "/level_%i",a_level);
  error = a_handle.setGroup(currentGroup + levelName);
  if (error != 0) return 1;

  HDF5HeaderData meta;
  meta.m_realvect["vec_dx"] = a_dx;
  meta.m_real["dt"] = a_dt;
  meta.m_real["time"] = a_time;
  meta.m_box["prob_domain"] = a_domain;
  meta.m_intvect["vec_ref_ratio"] = a_refRatios;

  error = meta.writeToFile(a_handle);
  if (error != 0) return 2;

  error = write(a_handle, a_data.boxLayout());
  if (error != 0) return 3;

  error = write(a_handle, a_data, "data", outputGhost, comps);
  if (error != 0) return 4;

  a_handle.setGroup(currentGroup);

  return 0;
}


template <class T>
int readLevel(HDF5Handle&   a_handle,
              const int&    a_level,
              LevelData<T>& a_data,
              Real& a_dx,
              Real& a_dt,
              Real& a_time,
              Box&  a_domain,
              int&  a_refRatio,
              const Interval&   a_comps,
              bool  setGhost)
{
  HDF5HeaderData header;
  header.readFromFile(a_handle);
  //unused
  // int nComp = header.m_int["num_components"];

  int error;
  char levelName[10];
  std::string currentGroup = a_handle.getGroup();
  sprintf(levelName, "/level_%i",a_level);
  error = a_handle.setGroup(currentGroup + levelName);
  if (error != 0) return 1;

  HDF5HeaderData meta;
  error = meta.readFromFile(a_handle);
  if (error != 0) return 2;
  a_dx       = meta.m_real["dx"];
  a_dt       = meta.m_real["dt"];
  a_time     = meta.m_real["time"];
  a_domain   = meta.m_box["prob_domain"];
  a_refRatio = meta.m_int["ref_ratio"];
  Vector<Box> boxes;
  error = read(a_handle, boxes);
  Vector<int> procIDs;
  LoadBalance(procIDs, boxes);

  DisjointBoxLayout layout(boxes, procIDs, a_domain);

  layout.close();
  if (error != 0) return 3;

  error = read(a_handle, a_data, "data", layout, a_comps, true);

  if (error != 0) return 4;

  a_handle.setGroup(currentGroup);

  return 0;
}

template <class T>
int readLevel(HDF5Handle&   a_handle,
              const int&    a_level,
              LevelData<T>& a_data,
              RealVect& a_dx,
              Real& a_dt,
              Real& a_time,
              Box&  a_domain,
              IntVect&  a_refRatio,
              const Interval&   a_comps,
              bool  setGhost)
{
  HDF5HeaderData header;
  header.readFromFile(a_handle);
  //unused
  // int nComp = header.m_int["num_components"];

  int error;
  char levelName[10];
  std::string currentGroup = a_handle.getGroup();
  sprintf(levelName, "/level_%i",a_level);
  error = a_handle.setGroup(currentGroup + levelName);
  if (error != 0) return 1;

  HDF5HeaderData meta;
  error = meta.readFromFile(a_handle);
  if (error != 0) return 2;
  a_dx       = meta.m_realvect["vec_dx"];
  a_dt       = meta.m_real["dt"];
  a_time     = meta.m_real["time"];
  a_domain   = meta.m_box["prob_domain"];
  a_refRatio = meta.m_intvect["vec_ref_ratio"];
  Vector<Box> boxes;
  error = read(a_handle, boxes);
  Vector<int> procIDs;
  LoadBalance(procIDs, boxes);

  DisjointBoxLayout layout(boxes, procIDs, a_domain);

  layout.close();
  if (error != 0) return 3;

  error = read(a_handle, a_data, "data", layout, a_comps, true);

  if (error != 0) return 4;

  a_handle.setGroup(currentGroup);

  return 0;
}

#include "NamespaceFooter.H"

#else // CH_USE_HDF5

// this is the only thing needed when HDF is not used
#define HOFFSET(S,M)    (offsetof(S,M))

#endif // CH_USE_HDF5 not defined
#endif  // CH_HDF5_H

---