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
// ============


template <class T>
int writeLevel(HDF5Handle& a_handle,
               const int& a_level,
               const LevelData<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);

// 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.


int write(HDF5Handle& a_handle,
          const BoxLayout& a_layout,
          const std::string& name = "boxes");


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


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());


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


int readBoxes(HDF5Handle& a_handle,
              Vector<Vector<Box> >& boxes);


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


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


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


class HDF5Handle
{
public:

  enum mode {CREATE, OPEN_RDONLY, OPEN_RDWR};



  HDF5Handle();


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

  ~HDF5Handle();



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


  bool isOpen() const;


  void close();



  void setGroupToLevel(int a_level);


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


  int pushGroup( const std::string& grp );


  int popGroup();


  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();

};


class HDF5HeaderData
{
public:


  int writeToFile(HDF5Handle& file) const;


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

  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)
{
  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);

  if(newForm && memcmp(a_name.c_str(), "M", 1)==0)
    {
      std::string currentGroup = a_handle.getGroup();
      a_handle.setGroup("/");
      HDF5HeaderData info;
      info.readFromFile(a_handle);
      //outputGhost = info.m_intvect["Ghost"];
      info.clear();
      a_handle.setGroup("/" + HDF5Handle::groups[a_name] + "CenteredComponents");
      info.readFromFile(a_handle);
      int nc = info.m_int["Num"+a_name];
      if(nc == 0){
        char cname[40];
        info.m_int["Num"+a_name] = comps.size();
        nc = comps.size();
        for(int i=0; i<nc; ++i){
          sprintf(cname, "Component%i", i);
          info.m_string[cname] = cname;
        }
        info.writeToFile(a_handle);
      }
      a_handle.setGroup(currentGroup);
      CH_assert(nc == comps.size());
    }

  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(memcmp(a_name.c_str(), "M", 1)==0){
          strcpy(dataname, "Mask");
        }else{
          sprintf(dataname, "%sRegular", a_name.c_str());
        }
      }else {
        sprintf(dataname, "%s:datatype=%i",a_name.c_str(), i);
      }
      dataspace[i]      = H5Screate_simple(1, flatdims, NULL);
      CH_assert(dataspace[i] >=0);
      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){
        sprintf(dataname, "%sOffsets",a_name.c_str());
      } else {
        sprintf(dataname, "%s:offsets=%i",a_name.c_str(), i);
      }
      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)
        {
          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);
        }
      H5Sclose(offsetspace);
      H5Dclose(offsetData);
    }

  // write BoxLayoutData attributes into Dataset[0]
  if(!newForm){
    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";
  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)
    {
      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
  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);
      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];
          if(count[0] > 0){
          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);
#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)
    {
      freeMT(buffers[i]);
      H5Sclose(dataspace[i]);
      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)
{
  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());
  for(DataIterator it = a_data.dataIterator(); it.ok(); ++it)
    {
      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)
        {
// Quick fix for Ted's tempest problem.  Brian will think about this.
//          if(a_comps.size() == 0){
            offset[0] = offsets[i][index];
            count[0] = offsets[i][index+1] - offset[0];
//          } else {
//            offset[0] = offsets[i][index] + box.numPts()*a_comps.begin();
//            count[0]  = a_comps.size() * box.numPts();
//          }
          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);
    }

 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 LevelData<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 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);

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

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

  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

---