TestHDF5.C

Go to the documentation of this file.

#include "Testing.H"
#include "Simulation.H"
#include "Stencil.H"
#include "PCPPCommUtil.H"
#include "PCPPIO.H"
#include "PCPPReport.H"
#include "PCPPIntervalUtils.H"
#include "EulerRHS.H"
#include "TestFixtures.H"
#include "EulerTestFixtures.H"
#include "PC2Util.H"
#include "Report.H"

using pcpp::comm::CheckResult;
typedef pcpp::ParallelGlobalType                   global_t;

void TestHDF5LegacyFileInfo(ix::test::results &parallelUnitResults,
                            pcpp::CommunicatorType &testComm)
{
  int myRank = testComm.Rank();
  pcpp::io::simfileinfo simFileInfo;  
  std::ostringstream messageStream;
  
  std::string testFileName2D("HDF5Examples/MultipleGrids2D.h5");
  std::string testFileName3D_1("HDF5Examples/SingleGrid3D.h5");
  std::string testFileName3D_2("HDF5Examples/MultipleGrids3D.h5");
  
  std::vector<size_t> expectedSizes1(2,61);
  std::vector<std::vector<size_t> > expectedSizes(2,expectedSizes1);
  expectedSizes[1][0] = 121;
  expectedSizes[1][1] = 40;

  std::bitset<pcpp::io::NUMFORMATBITS> expectedFormatBits;
  expectedFormatBits.reset();
  expectedFormatBits.set(pcpp::io::ISLEGACY);
  expectedFormatBits.set(pcpp::io::HASIBLANK);
  expectedFormatBits.set(pcpp::io::HASSTATE);
  expectedFormatBits.set(pcpp::io::HASTARGDATA);
  expectedFormatBits.set(pcpp::io::HASGRID);

  bool serialgetfileinfo2 = true;
  if(myRank == 0){
    if(pcpp::io::hdf5::FileInfo(testFileName2D,simFileInfo,messageStream)){
      std::cout << "pcpp::io::hdf5::FileInfo failed: " << messageStream.str() 
                << std::endl;
      testComm.SetExit(1);
    }
    pcpp::report::SimFileInfo(std::cout,simFileInfo);
    if(simFileInfo.formatBits != expectedFormatBits){
      serialgetfileinfo2 = false;
      testComm.SetExit(1);
    }
  }
  if(testComm.Check()){
    std::cout << "TestHDF5LegacyFileInfo::Failed!! Skipping tests." << std::endl;
    return;
  }

  testComm.StreamBroadCast(simFileInfo);

  bool parallelgetfileinfo2 = true;
  if(simFileInfo.numGrids != 2 || (simFileInfo.formatBits != expectedFormatBits)){
    std::cout << "TestHDF5LegacyFileInfo:Error: file info wrong." << std::endl;
    parallelgetfileinfo2 = false;
  } else {
    
    for(int iGrid = 0;iGrid < simFileInfo.numGrids;iGrid++){
      std::string &gridName(simFileInfo.gridNames[iGrid]);
      std::ostringstream expectedGridName;
      expectedGridName << "Group00" << iGrid+1;
      
      if(gridName != expectedGridName.str()){
        parallelgetfileinfo2 = false;
      }
      
      if(simFileInfo.gridNumDimensions[iGrid] != 2){
        parallelgetfileinfo2 = false;
      }
      
      std::vector<size_t> &gridSizes(simFileInfo.gridSizes[iGrid]);
      //      fixtures::IndexIntervalType &gridInterval(simFileInfo.gridSizes[iGrid]);
      if(gridSizes != expectedSizes[iGrid]){
        parallelgetfileinfo2 = false;
      }
      
      if(myRank == 0){
        if(!parallelgetfileinfo2)
          serialgetfileinfo2 = false;
      }
    }
    if(myRank == 0)
      pcpp::report::SimFileInfo(std::cout,simFileInfo);
  }
  
  serialgetfileinfo2   = CheckResult(serialgetfileinfo2,testComm);
  parallelgetfileinfo2 = CheckResult(parallelgetfileinfo2,testComm);
  
  parallelUnitResults.UpdateResult("HDF5:Serial2DFileInfo",serialgetfileinfo2);
  parallelUnitResults.UpdateResult("HDF5:Parallel2DFileInfo",parallelgetfileinfo2);
  
  std::vector<size_t> expectedSizes3D1(3,21);
  // std::vector<fixtures::IndexIntervalType> expectedIntervals3(2);
  // expectedIntervals3[0].InitSimple(expectedSizes3);
  //  expectedIntervals2[1].InitSimple(expectedSizes2);

  bool serialgetfileinfo3 = true;
  bool parallelgetfileinfo3 = true;
  if(myRank == 0){
    if(pcpp::io::hdf5::FileInfo(testFileName3D_1,simFileInfo,messageStream)){
      std::cout << "pcpp::io::hdf5::FileInfo failed: " << messageStream.str() 
                << std::endl;
      serialgetfileinfo3 = false;
      testComm.SetExit(1);
    }
    pcpp::report::SimFileInfo(std::cout,simFileInfo);
    if(simFileInfo.formatBits != expectedFormatBits){
      serialgetfileinfo2 = false;
      testComm.SetExit(1);
    }
  }
  if(testComm.Check())
    return;
  testComm.StreamBroadCast(simFileInfo);

  if(simFileInfo.numGrids != 1 || (simFileInfo.formatBits != expectedFormatBits)){
    parallelgetfileinfo3 = false;
  } else {

    std::string &gridName(simFileInfo.gridNames[0]);
    std::string expectedGridName("Group001");
      
    if(gridName != expectedGridName){
      parallelgetfileinfo3 = false;
    }
    int numgDim = 3;  
    if(simFileInfo.gridNumDimensions[0] != numgDim){
      parallelgetfileinfo3 = false;
    }
      
    //    fixtures::IndexIntervalType &gridInterval(simFileInfo.gridSizes[0]);
    for(int iDim = 0;iDim < numgDim;iDim++){
      if(simFileInfo.gridSizes[0][iDim] != expectedSizes3D1[iDim]){
        parallelgetfileinfo3 = false;
      }
    }
    
    if(myRank == 0){
      if(!parallelgetfileinfo3)
        serialgetfileinfo3 = false;
    }
  }

  serialgetfileinfo3   = CheckResult(serialgetfileinfo3,testComm);
  parallelgetfileinfo3 = CheckResult(parallelgetfileinfo3,testComm);

  parallelUnitResults.UpdateResult("HDF5:Serial3DFileInfo",serialgetfileinfo3);
  parallelUnitResults.UpdateResult("HDF5:Parallel3DFileInfo",parallelgetfileinfo3);
  

  std::vector<size_t> expectedSizes3D2(3,20);
  expectedSizes3D2[0] = 41;

  bool serialgetfileinfo32 = true;
  bool parallelgetfileinfo32 = true;
  if(pcpp::io::hdf5::FileInfo(testFileName3D_2,simFileInfo,testComm,messageStream))
    return;

  //   if(myRank == 0){
  //     if(pcpp::io::hdf5::FileInfo(testFileName3D_2,simFileInfo,messageStream)){
  //       std::cout << "pcpp::io::hdf5::FileInfo failed: " << messageStream.str() 
  //                 << std::endl;
  //       serialgetfileinfo32 = false;
  //       testComm.SetExit(1);
  //     }
  //     pcpp::report::SimFileInfo(std::cout,simFileInfo);
  //     if(simFileInfo.formatBits != expectedFormatBits){
  //       serialgetfileinfo2 = false;
  //       testComm.SetExit(1);
  //     }
  //   }
  //   if(testComm.Check())
  //     return;
  //   testComm.StreamBroadCast(simFileInfo);
  
  if(simFileInfo.numGrids != 2 || (simFileInfo.formatBits != expectedFormatBits)){
    parallelgetfileinfo32 = false;
  } else {

    for(int iGrid = 0;iGrid < simFileInfo.numGrids;iGrid++){
      std::string &gridName(simFileInfo.gridNames[iGrid]);
      std::ostringstream expectedGridName;
      expectedGridName << "Group00" << iGrid+1;
      
      if(gridName != expectedGridName.str()){
        parallelgetfileinfo32 = false;
      }
      
      if(simFileInfo.gridNumDimensions[iGrid] != 3){
        parallelgetfileinfo32 = false;
      }
      
      //      fixtures::IndexIntervalType &gridInterval(simFileInfo.gridSizes[iGrid]);
      if(simFileInfo.gridSizes[iGrid] != expectedSizes3D2){
        parallelgetfileinfo32 = false;
      }
      
      if(myRank == 0){
        if(!parallelgetfileinfo32)
          serialgetfileinfo32 = false;
      }
    }
  }
  
  serialgetfileinfo32   = CheckResult(serialgetfileinfo32,testComm);
  parallelgetfileinfo32 = CheckResult(parallelgetfileinfo32,testComm);
  
  parallelUnitResults.UpdateResult("HDF5:Multi3DFileInfo",serialgetfileinfo32);
  parallelUnitResults.UpdateResult("HDF5:ParMulti3DFileInfo",parallelgetfileinfo32);
  
};

void TestHDF5Read2DBlockStructuredGrid(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{

  std::ostringstream messageStream;
  std::string testFileName2D("HDF5Examples/MultipleGrids2D.h5");

  pcpp::CommunicatorType myCartComm;  
  pcpp::ParallelTopologyInfoType cartInfo;
  cartInfo.cartDecompDirections.resize(2,1);
  cartInfo.numDimensions = 2;
  pcpp::comm::SetupCartesianTopology(testComm,cartInfo,myCartComm,messageStream);

  pcpp::report::CartesianSetup(messageStream,cartInfo);
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);

  std::vector<int> &cartCoords(myCartComm.CartCoordinates());
  std::vector<int> &cartDims(myCartComm.CartDimensions());
  std::vector<int> cartNeighbors;
  myCartComm.CartNeighbors(cartNeighbors);


  int myRank = testComm.Rank();
  int nProc  = testComm.Size();

  fixtures::ConfigurationType testConfig;
  typedef simulation::grid::parallel_blockstructured grid_t;
  pcpp::io::hdf5::configuration hdf5Config;  

  pcpp::io::simfileinfo simFileInfo;    
  if(myRank == 0){
    if(pcpp::io::hdf5::FileInfo(testFileName2D,simFileInfo,messageStream)){
      std::cout << "pcpp::io::hdf5::FileInfo failed: " << messageStream.str() 
                << std::endl;
      testComm.SetExit(1);
    }
  }
  if(testComm.Check()){
    std::cout << "TestHDF5Read2DBSGrid:ERROR!! FileInfo failed. Skipping tests." 
              << std::endl;
    return;
  }
  testComm.StreamBroadCast(simFileInfo);
  

  hdf5Config.collectiveMode = 1;
  int numGrids = simFileInfo.numGrids;
  for(int iGrid = 0;iGrid < simFileInfo.numGrids;iGrid++){

    std::string &gridName(simFileInfo.gridNames[iGrid]);

    grid_t testGrid;
    testGrid.SetGridSizes(simFileInfo.gridSizes[iGrid]);
    pcpp::IndexIntervalType globalExtent;
    pcpp::IndexIntervalType &partitionExtent(testGrid.PartitionInterval());
    
    globalExtent.InitSimple(simFileInfo.gridSizes[iGrid]);
    messageStream << "TestHDF5Read2DBSGrid:Global size: ";
    globalExtent.PrettyPrint(messageStream);
    messageStream << std::endl;
    pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
    pcpp::io::RenewStream(messageStream);

    std::cout << "TestHDF5Read2DBSGrid:Partitioning " << std::endl;
    if(pcpp::util::PartitionCartesianExtent(globalExtent,cartDims,cartCoords,partitionExtent,messageStream)){
      std::cout << "TestHDF5Read2DBSGrid: Partitioning failed for grid "
                    << iGrid << ". Skipping tests." << std::endl;
      return;
    }
    
    messageStream << "Partition interval: ";
    partitionExtent.PrettyPrint(messageStream);
    messageStream << std::endl;
    
    pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
    pcpp::io::RenewStream(messageStream);
    
    testComm.Barrier();
    testGrid.Finalize();
    testGrid.AllocateCoordinateData();
    
    if(pcpp::io::hdf5::ParallelReadBlockStructuredGrid<grid_t>(gridName,"/",testGrid,hdf5Config,testFileName2D,
                                                               testComm,messageStream)){
      std::cout << "pcpp::io::hdf5::ParallelReadBlockStructuredGrid failed: "
                << messageStream.str() << std::endl;
      testComm.SetExit(1);
    } else {
      pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
      pcpp::io::RenewStream(messageStream);
    }
  }                                                               
  bool result = true;
  if(testComm.Check())
    result = false;

  parallelUnitResults.UpdateResult("HDF5:Basic:testran",result);
  
}

void TestPCPPHDF5Read(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{

  std::cout << "TestPCPPHDF5Read Start!" << std::endl;

  std::ostringstream messageStream;
  std::string testFileName2D("HDF5Examples/MultipleGrids2D.h5");
  pcpp::io::hdf5::base hdf5File;
  
  int myRank = testComm.Rank();
  

  bool serialOpen          = true;
  bool serialOpenGroup     = true;
  bool serialOpenGroupFail = false;
  bool globalExistence     = true;
  bool scopeExistence      = true;
  bool readAttribute       = true;
  bool attExists           = true;
  bool readAttributeVec    = true;
  bool readLegacyHeader    = true;
  bool legacyHeaderExists  = true;
  bool dataDims            = true;
  bool attDims             = true;

  if(myRank == 0){

    if(hdf5File.Open(testFileName2D,true)){
      serialOpen = false;
    }
    if(!hdf5File.Exists("Group001"))
      scopeExistence = false;
    if(!hdf5File.Exists("/Group001"))
      globalExistence = false;
    if(!hdf5File.Exists("Group002/X"))
      scopeExistence = false;
    if(!hdf5File.Exists("/Group002/X"))
      globalExistence = false;

    hdf5File.AttributeExists("Group001/HEADER");
    hdf5File.AttributeExists("/Group001/HEADER");
    hdf5File.AttributeExists("Group001/useIB");
    hdf5File.AttributeExists("Group001/numberOfAuxVars");
    hdf5File.AttributeExists("Group002/gridSize");

    std::vector<size_t> groupHeaderDims(hdf5File.AttributeDimensions("/Group001/HEADER"));
    if(groupHeaderDims.size() != 1)
      attDims = false;
    else if (groupHeaderDims[0] != 4)
      attDims = false;

    std::vector<size_t> coordinateDimensions1(hdf5File.DataDimensions("/Group001/X"));
    std::vector<size_t> coordinateDimensions2(hdf5File.DataDimensions("/Group002/X"));
    if(coordinateDimensions1.size() != 2 ||
       coordinateDimensions2.size() != 2){
      dataDims = false;
    } else {
      if(coordinateDimensions1[0] != 61 ||
         coordinateDimensions1[1] != 61){
        dataDims = false;
      }
      if(coordinateDimensions2[0] != 40 ||
         coordinateDimensions2[1] != 121){
        dataDims = false;
      }
    }

    //    std::cout << "### Before number of Grids attribute exists" << std::endl;
      

    int numberOfGrids = 0;
    if(!hdf5File.AttributeExists("numberOfGrids"))
      attExists = false;
    if(hdf5File.ReadAttribute("numberOfGrids",&numberOfGrids))
      readAttribute = false;
    if(numberOfGrids != 2)
      readAttribute = false;
   
    double time = -1.0;
    double headerSpace[] = {-1,-1,-1,-1};
    //if(!hdf5File.AttributeExists("/HEADER"))
    if(!hdf5File.AttributeExists("HEADER")){
      legacyHeaderExists = false;
    }

    if (!readLegacyHeader) {
      std::cout << "Failed to correctly read the legacy PlasComCM HEADER" << std::endl;
      std::cout << "This test is written to always fail. It is currently broken on Vulcan" << std::endl;
    }

    // This should pass on all platforms, but it doesn't for some reason on Vulcan
    // H5ARead fails to read the time from the HEADER attribute, causing
    // All of the HDF5 tests to fail.
    if(legacyHeaderExists){
      if(hdf5File.ReadAttribute("HEADER",headerSpace))
        readAttribute = false;
      time = headerSpace[0];
      if(time != 0.0)
        readAttribute = false;
    } else {
      readAttribute = false;
    }

    std::vector<double> gridHeader(4,-1);
    hdf5File.OpenGroup("Group001");
    if(!hdf5File.AttributeExists("HEADER"))
      attExists = false;
    if(attExists){
      if(hdf5File.ReadAttribute("HEADER",&gridHeader[0]))
        readAttribute = false;
      if(gridHeader[0] != 0 ||
         gridHeader[1] != 0 ||
         gridHeader[2] != 0 ||
         gridHeader[3] != 0){
        readAttribute = false;
      }
    }
    hdf5File.CloseGroup();

    // open non-existent group
    if(hdf5File.OpenGroup("nope")){
      serialOpenGroupFail = true;
    }

    if(!serialOpenGroupFail)
      serialOpenGroup = false;

    // test for root group existence
    // if(!hdf5File.Exists("/")){
    //   globalExistence = false;
    // }

    // test for non-existent group
    if(hdf5File.Exists("Configuration")){
      globalExistence = false;
    }

    if(hdf5File.OpenGroup("Group001"))
      serialOpenGroup = false;

    if(!hdf5File.Exists("X")){
      scopeExistence = false;
    }

    // Should not resolve Group002 in relative path
    if(hdf5File.Exists("Group002")){
      scopeExistence = false;
    }

    // *Should* resolve /Group002 absolute path
    if(!hdf5File.Exists("/Group002")){
      globalExistence = false;
    }
    
    // Same relative/absolute path distinction for attribute
    if(hdf5File.AttributeExists("numberOfGrids"))
      attExists = false;

    if(!hdf5File.AttributeExists("/numberOfGrids"))
      attExists = false;

    std::cout << "ReadAttribute using C-style array." << std::endl;
    std::vector<size_t> gridSize(3,0);
    if(hdf5File.ReadAttribute("/Group001/gridSize",&gridSize[0]))
      readAttribute = false;
    if(gridSize[0] != 61 ||
       gridSize[1] != 61)
      readAttribute = false;

    std::cout << "ReadAttribute using STL Vector." << std::endl;
    std::vector<size_t> testSizeVec;
    if(hdf5File.ReadAttribute("/Group001/gridSize",testSizeVec))
      readAttributeVec = false;
    if(testSizeVec.size() != 2){
      readAttributeVec = false;
    } else if(testSizeVec[0] != 61 ||
              testSizeVec[1] != 61){
      readAttributeVec = false;
    }
        
    std::cout << "Closing HDF5 file." << std::endl;
    hdf5File.Close();
  }
  
  
  bool parallelOpen = true;
  if(hdf5File.Open(testFileName2D,testComm,true)){
    testComm.SetExit(1);
  }
  if(testComm.Check())
    parallelOpen = false;
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);

  parallelUnitResults.UpdateResult("HDF5:Object:SerialOpen",serialOpen);
  parallelUnitResults.UpdateResult("HDF5:Object:SerialOpenGroup",serialOpenGroup);
  parallelUnitResults.UpdateResult("HDF5:Object:ParallelOpen",parallelOpen);
  parallelUnitResults.UpdateResult("HDF5:Object:GlobalExistence",globalExistence);
  parallelUnitResults.UpdateResult("HDF5:Object:ScopeExistence",scopeExistence);
  parallelUnitResults.UpdateResult("HDF5:Object:AttributeExists",attExists);
  parallelUnitResults.UpdateResult("HDF5:Object:DataDimensions",dataDims);
  parallelUnitResults.UpdateResult("HDF5:Object:AttributeDimensions",attDims);
  parallelUnitResults.UpdateResult("HDF5:Object:ReadAttribute",readAttribute);
  parallelUnitResults.UpdateResult("HDF5:Object:ReadAttVector",readAttributeVec);
  parallelUnitResults.UpdateResult("HDF5:Object:LegacyHeaderExists",legacyHeaderExists);
  parallelUnitResults.UpdateResult("HDF5:Object:ReadLegacyHeader",readLegacyHeader);

};


void TestPCPPHDF5Write(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{

  int myRank = testComm.Rank();

  bool groupCreate    = true;
  bool attCreate      = true;
  bool dataCreate     = true;
  bool stringAttRead  = true;
  bool readDataSetVec = true;
  bool dataDims       = true;
  bool configIO       = true;

  std::string testFunctionName("TestPCPPHDF5Write");
  pcpp::ConfigType testConfig;
  std::ostringstream configStream;

  if(pcpp::util::ProcessConfigurationFile("Config/Examples/misc/advect1d.config",testConfig,testComm,std::cout)){
    configIO = false;
    std::cout << testFunctionName 
              << ":Error:Failed to process configuration file Config/Examples/misc/advect1d.config." << std::endl;
  } else {
    configStream << testConfig;
  }

  if(myRank == 0){

    if(ix::sys::FILEEXISTS("tmptest.h5"))
      ix::sys::Remove("tmptest.h5");

    pcpp::io::hdf5::base testFile("tmptest.h5");
    testFile.CreateGroup("/Geometry");
    testFile.CreateGroup("/Domains");
    testFile.CreateGroup("/Geometry/Grids");
    testFile.CreateGroup("/Configuration");
    
    std::string configString(configStream.str());
    if(!configString.empty())
      testFile.CreateAttribute("/Configuration/Parameters",configString);

    int numGrids = 2;
    std::vector<size_t> attSize2(1,1);
    testFile.CreateAttribute("/Geometry/Grids/numGrids",attSize2,&numGrids);
    std::string gridNames("Grid001 Grid002");
    testFile.CreateAttribute("/Geometry/Grids/gridNames",gridNames);
    testFile.CreateGroup("/Geometry/Grids/Grid001");
    testFile.CreateGroup("/Geometry/Grids/Grid002");

    int numDim = 2;
    std::vector<size_t> gridSize(2,2);
    gridSize[0] = 3;
    std::vector<size_t> sizeAtt(1,2);
    testFile.CreateAttribute("/Geometry/Grids/Grid001/numDim",attSize2,&numDim);
    testFile.CreateAttribute("/Geometry/Grids/Grid002/numDim",attSize2,&numDim);
    testFile.CreateAttribute("/Geometry/Grids/Grid001/gridSize",sizeAtt,&gridSize[0]);
    testFile.CreateAttribute("/Geometry/Grids/Grid002/gridSize",sizeAtt,gridSize);

    int numGridPoints = 1;
    for(int iDim = 0;iDim < numDim;iDim++)
      numGridPoints *= gridSize[iDim];
    std::vector<double> grid1(numGridPoints,0);
    std::vector<double> grid2(numGridPoints,0);
    for(int iPoint = 0;iPoint < numGridPoints;iPoint++){
      grid1[iPoint] = iPoint;
      grid2[iPoint] = numGridPoints + iPoint;
    }

    testFile.CreateDataSet("/Geometry/Grids/Grid001/X",gridSize,grid1);
    testFile.CreateDataSet("/Geometry/Grids/Grid001/Y",gridSize,grid1);
    testFile.CreateDataSet("/Geometry/Grids/Grid002/X",gridSize,grid2);
    testFile.CreateDataSet("/Geometry/Grids/Grid002/Y",gridSize,grid1);

    testFile.Close();

    testFile.Open("tmptest.h5",true);

    if(!testFile.Exists("/Geometry")               ||
       !testFile.Exists("/Configuration")          ||
       !testFile.Exists("/Geometry/Grids")         ||
       !testFile.Exists("/Geometry/Grids/Grid001") ||
       !testFile.Exists("/Geometry/Grids/Grid002") ||
       !testFile.Exists("/Domains"))
      groupCreate = false;
    
    if(!testFile.AttributeExists("/Geometry/Grids/numGrids")         ||
       !testFile.AttributeExists("/Geometry/Grids/gridNames")        ||
       !testFile.AttributeExists("/Geometry/Grids/Grid001/numDim")   ||
       !testFile.AttributeExists("/Geometry/Grids/Grid002/numDim")   ||
       !testFile.AttributeExists("/Geometry/Grids/Grid001/gridSize") ||
       !testFile.AttributeExists("/Geometry/Grids/Grid002/gridSize"))
      attCreate = false;
    
    if(!testFile.Exists("/Geometry/Grids/Grid001/X") ||
       !testFile.Exists("/Geometry/Grids/Grid002/X") ||
       !testFile.Exists("/Geometry/Grids/Grid001/Y") ||
       !testFile.Exists("/Geometry/Grids/Grid002/Y"))
      dataCreate = false;

    std::string gridNamesAtt;
    testFile.ReadAttribute("/Geometry/Grids/gridNames",gridNamesAtt);
    if(gridNamesAtt != "Grid001 Grid002"){
      std::cout << "Read GridNames = (" << gridNamesAtt 
                << ")" << std::endl;
      stringAttRead = false;
    }
    
    pcpp::ConfigType testConfig2;
    std::string configString2;
    testFile.ReadAttribute("/Configuration/Parameters",configString2);
    std::istringstream configStream(configString2);
    configStream >> testConfig2;
    if(testConfig2 != testConfig){
      std::cout << "======= CONFIGURATION 1 ========" << std::endl
                << testConfig << std::endl
                << "======= CONFIGURATION 2 ========" << std::endl
                << testConfig2 << std::endl;
      configIO = false;
    }

    std::vector<size_t> grid1Size(testFile.DataDimensions("/Geometry/Grids/Grid001/X"));
    if(grid1Size.size() != 2)
      dataDims = false;
    else {
      if(grid1Size[0] != 3 ||
         grid1Size[1] != 2){
        dataDims = false;
        std::cout << "SIZE PROBLEM: " << grid1Size[0]
                  << "," << grid1Size[1] << std::endl;
      }        
    }
    std::vector<double> grid1XData;
    if(testFile.ReadDataSet("/Geometry/Grids/Grid001/X",grid1XData))
      std::cout << "ERROR!!!" << std::endl;
    if(grid1XData.size() != numGridPoints)
      readDataSetVec = false;
    for(int iPoint = 0;iPoint < numGridPoints;iPoint++)
      if(grid1XData[iPoint] != iPoint)
        readDataSetVec = false;
    testFile.Close();
  }
  
  parallelUnitResults.UpdateResult("HDF5:Object:CreateGroup",groupCreate);
  parallelUnitResults.UpdateResult("HDF5:Object:CreateAttribute",attCreate);
  parallelUnitResults.UpdateResult("HDF5:Object:CreateData",dataCreate);
  parallelUnitResults.UpdateResult("HDF5:Object:ReadStringAttribute",stringAttRead);
  parallelUnitResults.UpdateResult("HDF5:Object:ReadDataVector",readDataSetVec);
  parallelUnitResults.UpdateResult("HDF5:Object:DataDims",dataDims);
  parallelUnitResults.UpdateResult("HDF5:Object:ConfigIO",configIO);

};

void TestHDF5HyperSlab(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{
  int myRank = testComm.Rank();
  std::ostringstream messageStream;

  int numDim = 3;
  std::vector<size_t> dataSize(numDim,0);
  dataSize[0] = 8;
  dataSize[1] = 4;
  dataSize[2] = 2;

  // REVERSAL DUE TO HDF5 AUTOMATIC DATA ORDERING
  // and CONVENTION OF PlasCom2 
  // HDF5 assumes that data is arranged with Z being
  // the stride-1 direction. PlasCom2 stores data
  // in contiguous linear buffers with X being the
  // stride-1 direction.
  std::vector<size_t> rdataSize(dataSize);
  std::reverse(rdataSize.begin(),rdataSize.end());

  size_t numGlobalPoints = 1;

  if(myRank == 0){
    if(ix::sys::FILEEXISTS("tmptest2.h5"))
      ix::sys::Remove("tmptest2.h5");
    pcpp::io::hdf5::base testFile("tmptest2.h5");
    for(int iDim = 0;iDim < numDim;iDim++)
      numGlobalPoints *= dataSize[iDim];
    std::vector<double> myData(numGlobalPoints,0);
    size_t iPoint = 0;
    for(size_t iZ = 0;iZ < dataSize[2];iZ++){
      for(size_t iY = 0;iY < dataSize[1];iY++){
        for(size_t iX = 0;iX < dataSize[0];iX++){
          myData[iPoint++] = iX+iY*100+iZ*1000;
        }
      }
    }
    testFile.CreateDataSet("/myData",rdataSize,myData);
    testFile.Close();
  }
  
  testComm.Barrier();

  pcpp::CommunicatorType myCartComm;  
  pcpp::ParallelTopologyInfoType cartInfo;
  pcpp::comm::SetupCartesianTopology(testComm,cartInfo,myCartComm,messageStream);

  pcpp::report::CartesianSetup(messageStream,cartInfo);
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);

  std::vector<int> &cartCoords(myCartComm.CartCoordinates());
  std::vector<int> &cartDims(myCartComm.CartDimensions());
  std::vector<int> cartNeighbors;
  myCartComm.CartNeighbors(cartNeighbors);

  pcpp::IndexIntervalType globalExtent;
  globalExtent.InitSimple(dataSize);
  pcpp::IndexIntervalType partitionExtent;
  
  if(myRank == 0){
    std::cout << "Global interval: ";
    globalExtent.PrettyPrint(std::cout);
    std::cout << std::endl;
  }

  if(pcpp::util::PartitionCartesianExtent(globalExtent,cartDims,cartCoords,partitionExtent,messageStream)){
    std::cout << "TestHDF5HyperSlab:ERROR: Partitioning failed." << std::endl;
    return;
  }
  
  messageStream << "Partition interval: ";
  partitionExtent.PrettyPrint(messageStream);
  messageStream << std::endl;
  
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);
  
  int numLocalPoints = partitionExtent.NNodes();
  std::vector<size_t> localStart(partitionExtent.Starts());
  std::vector<size_t> rlocalStart(localStart);
  std::reverse(rlocalStart.begin(),rlocalStart.end());
  std::vector<size_t> localSize(partitionExtent.Sizes());
  std::vector<size_t> rlocalSize(localSize);
  std::reverse(rlocalSize.begin(),rlocalSize.end());
  std::vector<double> myData;

  pcpp::io::hdf5::base testFile("tmptest2.h5",testComm);
  testFile.ReadHyperSlab("/myData",rlocalStart,rlocalSize,myData);

  bool readHyperSlab = true;
  if(myData.size() != numLocalPoints){
    readHyperSlab = false;
  } else {
    size_t iLocal = 0;
    for(size_t iZ = partitionExtent[2].first;iZ <= partitionExtent[2].second;iZ++){
      for(size_t iY = partitionExtent[1].first;iY <= partitionExtent[1].second;iY++){
        for(size_t iX = partitionExtent[0].first;iX <= partitionExtent[0].second;iX++){
          double expectedValue = iX + 100*iY + 1000*iZ;
          if(myData[iLocal] != expectedValue){
            messageStream << myRank << "(local,[ix,iy,iz]) = (actual,expected) : ("
                          << iLocal << ",[" << iX << "," << iY << "," << iZ << "]) = ("
                          << myData[iLocal] << "," << expectedValue << ")" << std::endl;
            readHyperSlab = false;
          }
          iLocal++;
        }
      }
    }
  }
  testFile.Close();
  bool writeHyperSlab = false;
  if(readHyperSlab){
    writeHyperSlab = true;
    if(myRank == 0){
      ix::sys::Remove("tmptest2.h5");
    }
    testComm.Barrier();
    testFile.Create("tmptest2.h5",testComm);
    for(size_t iLocal = 0;iLocal < numLocalPoints;iLocal++){
      myData[iLocal] = -1.0*myData[iLocal];
    }
    testFile.WriteHyperSlab("/myData",rdataSize,rlocalStart,rlocalSize,myData);
    testFile.Close();
    testFile.Open("tmptest2.h5",testComm);
    myData.resize(0);
    testFile.ReadHyperSlab("/myData",rlocalStart,rlocalSize,myData);
    testFile.Close();
    size_t iLocal = 0;
    for(size_t iZ = partitionExtent[2].first;iZ <= partitionExtent[2].second;iZ++){
      for(size_t iY = partitionExtent[1].first;iY <= partitionExtent[1].second;iY++){
        for(size_t iX = partitionExtent[0].first;iX <= partitionExtent[0].second;iX++){
          double expectedValue = -1.0*(iX + 100*iY + 1000*iZ);
          if(myData[iLocal] != expectedValue){
            messageStream << myRank << "(local,[ix,iy,iz]) = (actual,expected) : ("
                          << iLocal << ",[" << iX << "," << iY << "," << iZ << "]) = ("
                          << myData[iLocal] << "," << expectedValue << ")" << std::endl;
            writeHyperSlab = false;
          }
          iLocal++;
        }
      }
    }
    ix::sys::Remove("tmptest2.h5");
  }
  
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);

  parallelUnitResults.UpdateResult("HDF5:Object:ReadHyperSlab",readHyperSlab);
  parallelUnitResults.UpdateResult("HDF5:Object:WriteHyperSlab",writeHyperSlab);

  if(myRank == 0){
    if(ix::sys::FILEEXISTS("tmptest2.h5"))
      ix::sys::Remove("tmptest2.h5");
  }
};                                                               


void TestHDF5GridIO(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{

  typedef simulation::grid::parallel_blockstructured         pbsgrid_t;
  typedef simulation::state::base                            state_t;

  int myRank  = testComm.Rank();
  int numProc = testComm.Size();

  std::ostringstream messageStream;

  pbsgrid_t   euler3dGrid;
  state_t  simulationState;

  double   t = 0.0;
  double   gamma = 1.4;
  double   sigma = 0.0;
  double   cfl = 1.0;

  // references to some grid-specific data structures
  pcpp::ParallelTopologyInfoType cartInfo;  
  std::vector<int>        &decompDirs(cartInfo.cartDecompDirections);
  std::vector<int>        &periodicDirs(cartInfo.isPeriodic);
  std::vector<double>     &realGridExtent(euler3dGrid.PhysicalExtent());
  std::vector<size_t>     &numGlobalNodes(euler3dGrid.GridSizes());
  //  pcpp::IndexIntervalType &globalExtent(euler3dGrid.);

  realGridExtent.resize(6,0.0);
  realGridExtent[1] = 1.0;
  realGridExtent[3] = 1.0;
  realGridExtent[5] = 1.0;

  numGlobalNodes.resize(3);
  numGlobalNodes[0] = 16;
  numGlobalNodes[1] = 8;
  numGlobalNodes[2] = 4;

  decompDirs.resize(3,1);
  //  decompDirs[0] = 1;

  periodicDirs.resize(3,1);
  euler3dGrid.ParallelSetup(testComm,cartInfo,3,messageStream);
 
  // Grab the local coordinates and global dimension of the Cartesian topo
  pcpp::CommunicatorType  &gridComm(euler3dGrid.Communicator());
  std::vector<int>        &cartCoords(gridComm.CartCoordinates());
  std::vector<int>        &cartDims(gridComm.CartDimensions());
  pcpp::IndexIntervalType &partitionExtent(euler3dGrid.PartitionInterval());

  // Grab the +/- neighbors in each Cartesian direction
  std::vector<int> cartNeighbors;
  gridComm.CartNeighbors(cartNeighbors);

  // Generate a report of the Cartesian setup and put it to stdout
  pcpp::report::CartesianSetup(messageStream,cartInfo);
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);
  
  //    dx[0] = (realGridExtent[1] - realGridExtent[0])/static_cast<double>(numGlobalNodes[0]-1);

  // Basic setup for the grid and record some local information
  std::vector<double> dX(3,0.0);
  dX[0] = (realGridExtent[1] - realGridExtent[0])/(static_cast<double>(numGlobalNodes[0]-1));
  dX[1] = (realGridExtent[3] - realGridExtent[2])/(static_cast<double>(numGlobalNodes[1]-1));
  dX[2] = (realGridExtent[5] - realGridExtent[4])/(static_cast<double>(numGlobalNodes[2]-1));
  euler3dGrid.SetGridSpacings(dX);
  
  pcpp::IndexIntervalType globalExtent;
  globalExtent.InitSimple(numGlobalNodes); 

  size_t totalGlobalNodes = numGlobalNodes[0]*numGlobalNodes[1]*numGlobalNodes[2];
  size_t totalGlobalCells = (numGlobalNodes[0]-1)*(numGlobalNodes[1]-1)*(numGlobalNodes[2]-1);

  double dx = dX[0];
  double dy = dX[1];
  double dz = dX[2];
  double dt = cfl*dx;
  dt = std::min(dt,cfl*dy);
  dt = std::min(dt,cfl*dz);
  //  eulerDomain.numGrids = 1;
  
  

  messageStream << "Global size: ";
  globalExtent.PrettyPrint(messageStream);
  messageStream << std::endl
                << "Partition interval: ";
  partitionExtent.PrettyPrint(messageStream);
  messageStream << std::endl;
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);

  size_t numX    = partitionExtent.NNodes(0);
  size_t numY    = partitionExtent.NNodes(1);
  size_t numZ    = partitionExtent.NNodes(2);
  size_t nNodes  = partitionExtent.NNodes();
  size_t nCells  = (numX-1)*(numY-1)*(numZ-1);
  size_t xyPlane = numX*numY;

  const std::vector<size_t>     &numPartitionNodes(euler3dGrid.PartitionSizes());
  // =====================

  // Hand the grid to the domain
  //  eulerDomain.Grids().resize(1,&euler3dGrid);

  // This input state is just so that
  // the domain object can read the meta
  // data and sizes.  The domain object 
  // allocates its own version of this state.
  simulationState.AddField("rho",'n',1,8,"mass");
  simulationState.AddField("rhoV",'n',3,8,"momentum");
  simulationState.AddField("rhoE",'n',1,8,"energy");
  simulationState.AddField("pressure",'n',1,8,"pressure");
  simulationState.AddField("temperature",'n',1,8,"temperature");
  simulationState.AddField("simTime",'s',1,8,"s");
  simulationState.AddField("inputDT",'s',1,8,"s");
  simulationState.AddField("inputCFL",'s',1,8,"");
  simulationState.AddField("gamma",'s',1,8,"");
  simulationState.AddField("refSigma",'s',1,8,"");


  // Specify which fields will advance (i.e. other than time)
  std::vector<int> stateFieldIndices(3,0);
  stateFieldIndices[1] = 1;
  stateFieldIndices[2] = 2;
  simulationState.SetStateFields(stateFieldIndices);
  simulationState.Create(nNodes,0);
  


  // === Set up the initial state ===
  size_t iNode = 0;
  std::vector<double> rhoVector;
  std::vector<double> rhoVVector;
  std::vector<double> rhoEVector;
  std::vector<double> tVector;
  std::vector<double> pVector;
  rhoVector.reserve(nNodes);
  rhoVVector.reserve(3*nNodes);
  rhoEVector.reserve(nNodes);
  tVector.reserve(nNodes);
  pVector.reserve(nNodes);
  double *rhoBuffer = &rhoVector[0];
  double *rhoVBuffer = &rhoVVector[0];
  double *rhoEBuffer = &rhoEVector[0];
  double *temperatureBuffer = &tVector[0];
  double *pressureBuffer    = &pVector[0];

  // Link the domain's state storage with local storage
  simulationState.SetFieldBuffer("rho",rhoBuffer);
  simulationState.SetFieldBuffer("rhoV",rhoVBuffer);
  simulationState.SetFieldBuffer("rhoE",rhoEBuffer);
  simulationState.SetFieldBuffer("temperature",temperatureBuffer);
  simulationState.SetFieldBuffer("pressure",pressureBuffer);
  simulationState.SetFieldBuffer("simTime",&t);
  simulationState.SetFieldBuffer("inputDT",&dt);
  simulationState.SetFieldBuffer("gamma",&gamma);
  simulationState.SetFieldBuffer("refSigma",&sigma);
  simulationState.SetFieldBuffer("inputCFL",&cfl);
  
  // Report the problem setup to stdout
  messageStream << "Problem setup:    " << std::endl
                << "-------------------------------"  << std::endl; 
  if(myRank == 0){
    messageStream << "Sigma:            " << sigma      << std::endl
                  << "CFL:              " << cfl        << std::endl
                  << "DX:               " << dx         << std::endl
                  << "DY:               " << dy         << std::endl
                  << "DZ:               " << dz         << std::endl
                  << "DT:               " << dt         << std::endl
                  << "Domain:           [" << realGridExtent[0] << "," 
                  << realGridExtent[1] << "] x [" << realGridExtent[2] 
                  << "," << realGridExtent[3] << "] x [" << realGridExtent[4] 
                  << "," << realGridExtent[5] << "]" << std::endl
                  << "Grid Points:      (" << numGlobalNodes[0] 
                  << "," << numGlobalNodes[1] << "," << numGlobalNodes[2] 
                  << ")" << std::endl;
  }
  messageStream << "Partition Interval: ";
  euler3dGrid.PartitionInterval().PrettyPrint(messageStream);
  messageStream << std::endl
                << "-------------------------------"  << std::endl << std::endl;
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);
  
  if(ix::sys::FILEEXISTS("gridtest.h5"))
    ix::sys::Remove("gridtest.h5");
  pcpp::io::hdf5::base gridFile("gridtest.h5",gridComm);
  gridFile.CreateGroup("/Global");
  gridFile.CreateGroup("/PlasCom2");
  gridFile.CreateGroup("/PlasCom2/Configuration");
  gridFile.CreateGroup("/PlasCom2/Geometry");
  gridFile.CreateGroup("/PlasCom2/Geometry/Grids");
  
  bool writeGridTest = true;
  int numGrids = 2;
  std::vector<size_t> attSize2(1,1);
  int numGeometries = 1;
  std::string geometryNames("Grids");
  gridFile.CreateAttribute("/PlasCom2/Geometry/numGeometries",attSize2,&numGeometries);
  gridFile.CreateAttribute("/PlasCom2/Geometry/geometryNames",geometryNames);
  gridFile.CreateAttribute("/PlasCom2/Geometry/Grids/numGrids",attSize2,&numGrids);
  std::string gridNames("Grid001 Grid002");
  gridFile.CreateAttribute("/PlasCom2/Geometry/Grids/gridNames",gridNames);
  std::vector<std::string> expectedGridNames(2,"Grids:Grid001");
  expectedGridNames[1] = "Grids:Grid002";
  
  pcpp::IndexIntervalType preDimension;
  preDimension.InitSimple(numGlobalNodes);
  std::cout << "Pre-Generation Dimension: ";
  preDimension.PrettyPrint(std::cout);
  std::cout << std::endl;

  //  euler3dGrid.InitDone();

  //  euler3dGrid.GenerateCoordinateData();
  euler3dGrid.GenerateCoordinates(std::cout);
  std::vector<double> expectedCoordinates(euler3dGrid.CoordinateData());

  // std::cout << "EXPECTED COORDS: " << std::endl;
  // pcpp::io::DumpContents(std::cout,expectedCoordinates,",");
  // std::cout << std::endl;
  
  preDimension.InitSimple(numGlobalNodes);
  std::cout << "Post-Generation Dimension: ";
  preDimension.PrettyPrint(std::cout);
  std::cout << std::endl;

  pcpp::io::hdf5::WriteGrid(euler3dGrid,"Grid001","/PlasCom2/Geometry/Grids",gridFile);
  pcpp::io::hdf5::WriteGrid(euler3dGrid,"Grid002","/PlasCom2/Geometry/Grids",gridFile);

  gridFile.Close();

  //  return;

  pcpp::io::simfileinfo simFileInfo;  
  if(myRank == 0){
    if(pcpp::io::hdf5::FileInfo("gridtest.h5",simFileInfo,messageStream)){
      std::cout << "pcpp::io::hdf5::FileInfo failed: " << messageStream.str() 
                << std::endl;
      testComm.SetExit(1);
    }
    pcpp::report::SimFileInfo(messageStream,simFileInfo);
  }
  if(testComm.Check()){
    std::cout << "TestHDF5GridIO:ERROR!! FileInfo failed. Skipping tests." << std::endl;
    return;
  }
  
  if(myRank == 0){
    std::ostringstream infoStreamOut;
    infoStreamOut << simFileInfo;
    messageStream << "SIMFILEINFO PACKED:" << std::endl
                  << "(" << infoStreamOut.str() << ")" << std::endl;
    std::istringstream infoStreamIn(infoStreamOut.str());
    pcpp::io::simfileinfo info2;
    infoStreamIn >> info2;
    messageStream << "SIMFILEINFO UnPACKED:" << std::endl
                  << "(" << info2 << ")" << std::endl;
  } 
  
  
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);

  testComm.StreamBroadCast(simFileInfo);

  if(myRank != 0){
    messageStream << "REmote Unpacked:" << std::endl
                  << "(" << simFileInfo << ")" << std::endl;
  }

  pcpp::report::SimFileInfo(messageStream,simFileInfo);
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);
  
  pcpp::IndexIntervalType expectedDimension;
  expectedDimension.InitSimple(numGlobalNodes);
  messageStream << "Input Grid Dimension: ";
  expectedDimension.PrettyPrint(std::cout);
  messageStream << std::endl;
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);

  if(simFileInfo.numGrids != 2){
    messageStream << "Unexpected number of grids: " << simFileInfo.numGrids 
                  << " != 2" << std::endl;
    writeGridTest = false;
  }
  if(simFileInfo.gridNames != expectedGridNames){
    messageStream << "Unpexected grid names: (";
    pcpp::io::DumpContents(messageStream,simFileInfo.gridNames,",");
    messageStream << ") != (";
    pcpp::io::DumpContents(messageStream,expectedGridNames,",");
    messageStream << ")" << std::endl;
    writeGridTest = false;
  }
  for(int iGrid = 0;iGrid < simFileInfo.numGrids;iGrid++){
    if(simFileInfo.gridNumDimensions[iGrid] != 3){
      messageStream << "Unexpected grid numDim for grid " << iGrid 
                    << ": " << simFileInfo.gridNumDimensions[iGrid] 
                    << " != 3" << std::endl;
      writeGridTest = false;
    }
    if(simFileInfo.gridSizes[iGrid] != numGlobalNodes){
      messageStream << "Unpexected grid size for grid " << iGrid 
                    << ": (";
      pcpp::io::DumpContents(messageStream,simFileInfo.gridSizes[iGrid],",");
      messageStream << ") != (";
      pcpp::io::DumpContents(messageStream,numGlobalNodes,",");
      messageStream << ")" << std::endl;
      writeGridTest = false;
    }
  }
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);

  writeGridTest = CheckResult(writeGridTest,testComm);
  parallelUnitResults.UpdateResult("HDF5:Interface:WriteGrid",writeGridTest);

  // Re-open read-only
  

  bool readGridTest = true;

  // reset coordinate data to zero 
  std::vector<double> &modifyCoordinateData(euler3dGrid.CoordinateData());
  std::vector<double>::iterator cdIt = modifyCoordinateData.begin();
  while(cdIt != modifyCoordinateData.end())
    *cdIt++ = 0;

  gridFile.Open("gridtest.h5",gridComm,true);
  if(pcpp::io::hdf5::ReadGrid(euler3dGrid,"/PlasCom2/Geometry/Grids/Grid001",gridFile)){
    std::cout << "Reading grid from gridtest.h5 failed!" << std::endl;
    readGridTest = false;
  }
  gridFile.Close();

  std::vector<double> &readCoordinates(euler3dGrid.CoordinateData());
  if(readCoordinates != expectedCoordinates){
    readGridTest = false;
    messageStream << "READ UNEXPECTED COORDS: " << std::endl;
    pcpp::io::DumpContents(messageStream,readCoordinates,",");
    messageStream << std::endl
                  << "--------WHICH IS NOT EQUAL TO EXPECTED:";
    pcpp::io::DumpContents(messageStream,expectedCoordinates,",");
    messageStream << std::endl;
  }
  pcpp::io::Everyone(messageStream.str(),std::cout,testComm);
  pcpp::io::RenewStream(messageStream);

  readGridTest = CheckResult(readGridTest,testComm);
  parallelUnitResults.UpdateResult("HDF5:Interface:ReadGrid",readGridTest);

  if(myRank == 0){
    if(ix::sys::FILEEXISTS("gridtest.h5"))
      ix::sys::Remove("gridtest.h5");
  }

};                                                               


void TestHDF5PBSGridIO(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{  
  std::ostringstream messageStream;

  fixtures::ConfigurationType testConfig;

  typedef simulation::grid::parallel_blockstructured pbsgrid_t;
  typedef plascom2::operators::sbp::base             operator_t;
  typedef simulation::grid::halo                     halo_t;
  typedef simulation::state::base                    state_t;
  typedef pcpp::CommunicatorType                     comm_t;
  typedef pcpp::IndexIntervalType                    interval_t;

  
  global_t myGlobal;
  myGlobal.Init("TestPBSIO",testComm);
  myGlobal.SetVerbLevel(3);
  myGlobal.Profiling(true);

  int numProc = testComm.Size();
  int myRank  = testComm.Rank();

  //  grid_t     testGrid;
  pbsgrid_t  pbsGrid;
  //  halo_t     testHalo;
  //  comm_t    &gridComm(testGrid.Communicator());
  int numDim = 3;
  std::vector<size_t> gridSize(numDim,0);
  std::vector<double> dX(numDim,0);

  size_t numGlobalNodes = 1;
  size_t numGlobalCells = 1;
  for(int iDim = 0;iDim < numDim;iDim++){
    // 100x50x25
    gridSize[iDim] = static_cast<size_t>(std::pow(2.0,(numDim-iDim-1)%numDim)*25);
    dX[iDim] = 1.0/static_cast<double>(gridSize[iDim]-1);
    numGlobalNodes *= gridSize[iDim];
    numGlobalCells *= (gridSize[iDim]-1);
  }
  // gridSize[0] = 4;
  // gridSize[1] = 4;
  // gridSize[2] = 2;
  // numGlobalNodes = 32;
  // numGlobalCells = 9;
  pbsGrid.SetGridSizes(gridSize);
  pcpp::ParallelTopologyInfoType pbsCartInfo;
  pbsCartInfo.numDimensions = numDim;
  pbsCartInfo.isPeriodic.resize(numDim,1);

  operator_t sbpOperator;
  plascom2::operators::sbp::Initialize(sbpOperator,4);
  int boundaryDepth = (sbpOperator.numStencils-1)/2;

  pbsGrid.ParallelSetup(testComm,pbsCartInfo,boundaryDepth,messageStream);

  std::vector<int> haloDepths(2*numDim,boundaryDepth);
  std::vector<int> boundaryDepths(2*numDim,boundaryDepth);

  comm_t    &pbsGridComm(pbsGrid.Communicator());
  // Grab the local coordinates and global dimension of the Cartesian topo
  std::vector<int> &cartCoords(pbsGridComm.CartCoordinates());
  std::vector<int> &cartDims(pbsGridComm.CartDimensions());
  // Grab the +/- neighbors in each Cartesian direction
  std::vector<int> cartNeighbors;
  pbsGridComm.CartNeighbors(cartNeighbors);

  pcpp::ParallelTopologyInfoType cartInfo;
  cartInfo.numDimensions = numDim;
  
  // Generate a report of the Cartesian setup and put it to stdout
  pcpp::report::CartesianSetup(messageStream,pbsCartInfo);
  myGlobal.StdOut(messageStream.str(),2);
  pcpp::io::RenewStream(messageStream);

  interval_t globalInterval;
  globalInterval.InitSimple(gridSize);

  //  interval_t &partitionInterval(pbsGrid.PartitionInterval());
  interval_t &pbsPartInterval(pbsGrid.PartitionInterval());
  const std::vector<size_t> &bufferSizes(pbsGrid.BufferSizes());
  const pcpp::IndexIntervalType &localPartitionInterval(pbsGrid.PartitionBufferInterval());  

  size_t numPointsBuffer = pbsGrid.BufferSize();

  pcpp::IndexIntervalType bufferInterval;
  bufferInterval.InitSimple(bufferSizes);


  // Create simple message buffers - assuming the grid coordinates will
  // be used to test 
  int numComponents = numDim;

  pbsGrid.AllocateCoordinateData();
  std::vector<double> &gridCoordinates(pbsGrid.CoordinateData());

  size_t partStartX   = pbsPartInterval[0].first;
  size_t partEndX     = pbsPartInterval[0].second;
  size_t partStartY   = pbsPartInterval[1].first;
  size_t partEndY     = pbsPartInterval[1].second;
  size_t partStartZ   = pbsPartInterval[2].first;
  size_t partEndZ     = pbsPartInterval[2].second;
  size_t bufferStartX = localPartitionInterval[0].first;
  size_t bufferEndX   = localPartitionInterval[0].second;
  size_t bufferStartY = localPartitionInterval[1].first;
  size_t bufferEndY   = localPartitionInterval[1].second;
  size_t bufferStartZ = localPartitionInterval[2].first;
  size_t bufferEndZ   = localPartitionInterval[2].second;
  size_t buffPlane    = bufferSizes[0]*bufferSizes[1];
  size_t buffRow      = bufferSizes[0];
  size_t iPartZ,iPartY,iPartX,iBufX, iBufY, iBufZ;
  for(iPartZ = partStartZ, iBufZ = bufferStartZ;
      iPartZ <= partEndZ;iPartZ++,iBufZ++){
    size_t bufZ = iBufZ*buffPlane;
    for(iPartY = partStartY, iBufY = bufferStartY;
        iPartY <= partEndY;iPartY++,iBufY++){
      size_t bufYZ = bufZ + iBufY*buffRow;
      for(iPartX = partStartX, iBufX = bufferStartX;
          iPartX <= partEndX;iPartX++,iBufX++){
        size_t bufIndex = bufYZ + iBufX;
        gridCoordinates[bufIndex]                   = iPartX;
        gridCoordinates[bufIndex+numPointsBuffer]   = iPartY;
        gridCoordinates[bufIndex+2*numPointsBuffer] = iPartZ;
      }
    }
  }
  pbsGrid.ExchangeCoordinates();
  //  testHalo.PackMessageBuffers(xyzMessageId,0,3,&gridCoordinates[0]);

  // Send Test w/Recv Buffer check
  //   testHalo.SendMessage(0,cartNeighbors,pbsGridComm);
  //   testHalo.ReceiveMessage(0,cartNeighbors,pbsGridComm);
  //   testHalo.UnPackMessageBuffers(0,0,3,&gridCoordinates[0]);
  
  if(myRank == 0){
    if(ix::sys::FILEEXISTS("pbsgridtest.h5"))
      ix::sys::Remove("pbsgridtest.h5");
  }
  testComm.Barrier();

  pcpp::io::hdf5::base gridFile("pbsgridtest.h5",pbsGridComm);

  gridFile.CreateGroup("/Global");
  gridFile.CreateGroup("/PlasCom2");
  gridFile.CreateGroup("/PlasCom2/Configuration");
  gridFile.CreateGroup("/PlasCom2/Geometry");
  gridFile.CreateGroup("/PlasCom2/Geometry/Grids");
  
  bool writeGridTest = true;
  std::vector<size_t> attSize2(1,1);
  
  int numGeometries = 1;
  std::string geometryNames("Grids");
  gridFile.CreateAttribute("/PlasCom2/Geometry/numGeometries",attSize2,&numGeometries);
  gridFile.CreateAttribute("/PlasCom2/Geometry/geometryNames",geometryNames);
  
  int numGrids = 1;
  gridFile.CreateAttribute("/PlasCom2/Geometry/Grids/numGrids",attSize2,&numGrids);
  std::string gridNames("pbsgrid");
  gridFile.CreateAttribute("/PlasCom2/Geometry/Grids/gridNames",gridNames);
  std::string expectedGridNames("Grids:pbsgrid");
  
  std::vector<double> expectedCoordinates(pbsGrid.CoordinateData());
  
  pcpp::io::hdf5::WriteGrid(pbsGrid,"pbsgrid","/PlasCom2/Geometry/Grids",gridFile);
  
  gridFile.Close();
  
  pcpp::io::simfileinfo simFileInfo;  
  if(myRank == 0){
    if(pcpp::io::hdf5::FileInfo("pbsgridtest.h5",simFileInfo,messageStream)){
      std::cout << "pcpp::io::hdf5::FileInfo failed: " << messageStream.str() 
                << std::endl;
      testComm.SetExit(1);
    }
    pcpp::report::SimFileInfo(std::cout,simFileInfo);
  }
  if(testComm.Check()){
    std::cout << "TestPBSGridIO:ERROR!! FileInfo failed. Skipping tests." << std::endl;
    return;
  }
  testComm.StreamBroadCast(simFileInfo);
  
  pcpp::IndexIntervalType expectedDimension;
  expectedDimension.InitSimple(gridSize);
  std::cout << "Input Grid Dimension: ";
  expectedDimension.PrettyPrint(std::cout);
  std::cout << std::endl;
  
  
  if(simFileInfo.numGrids != 1)
    writeGridTest = false;
  if(simFileInfo.gridNames[0] != expectedGridNames)
    writeGridTest = false;
  for(int iGrid = 0;iGrid < simFileInfo.numGrids;iGrid++){
    if(simFileInfo.gridNumDimensions[iGrid] != 3)
      writeGridTest = false;
    if(simFileInfo.gridSizes[iGrid] != gridSize)
      writeGridTest = false;
  }
  
  writeGridTest = CheckResult(writeGridTest,testComm);
  parallelUnitResults.UpdateResult("HDF5:Interface:WritePBSGrid",writeGridTest);
  
  
  bool readGridTest = true;
  
  // reset coordinate data to zero 
  std::vector<double> &modifyCoordinateData(pbsGrid.CoordinateData());
  std::vector<double>::iterator cdIt = modifyCoordinateData.begin();
  while(cdIt != modifyCoordinateData.end())
    *cdIt++ = 0;
  
  // Re-open read-only
  gridFile.Open("pbsgridtest.h5",pbsGridComm,true);
  int returnCode = pcpp::io::hdf5::ReadGrid(pbsGrid,"/PlasCom2/Geometry/Grids/pbsgrid",gridFile);
  if(returnCode > 0){
    std::cout << "Grid Reading failed for pbsgrid(" << returnCode 
              << ")." << std::endl;
    readGridTest = false;
  }
  gridFile.Close();
  pbsGrid.ExchangeCoordinates();
  
  // Send Test w/Recv Buffer check
  //   testHalo.PackMessageBuffers(xyzMessageId,0,3,&gridCoordinates[0]);
  //   testHalo.SendMessage(xyzMessageId,cartNeighbors,pbsGridComm);
  //   testHalo.ReceiveMessage(xyzMessageId,cartNeighbors,pbsGridComm);
  //   testHalo.UnPackMessageBuffers(xyzMessageId,0,3,&gridCoordinates[0]);
  
  //  std::vector<double> &readCoordinates(euler3dGrid.CoordinateData());
  if(gridCoordinates != expectedCoordinates)
    readGridTest = false;
  
  readGridTest = CheckResult(readGridTest,testComm);
  parallelUnitResults.UpdateResult("HDF5:Interface:ReadPBSGrid",readGridTest);

  state_t simulationState;
  state_t testState;
  // This input state is just so that
  // the domain object can read the meta
  // data and sizes.  The domain object 
  // allocates its own version of this state.
  simulationState.AddField("meshScalar",'m',1,8,"unit2");
  simulationState.AddField("meshVector0",'m',3,8,"");
  simulationState.AddField("meshVector1",'m',3,8,"unit1");
  simulationState.AddField("nodeScalar",'n',1,8,"unit2");
  simulationState.AddField("nodeVector0",'n',3,8,"");
  simulationState.AddField("nodeVector1",'n',3,8,"unit1");
  simulationState.AddField("nodeFlag",'n',1,4,"");
  simulationState.AddField("simVector0",'s',10,8,"");
  simulationState.AddField("simVector1",'s',2,8,"unit1");
  simulationState.AddField("simScalar",'s',1,8,"unit2");
  simulationState.AddField("simBytes",'s',10,1,"bytes");
  simulationState.AddField("simInts",'s',5,4,"ints");

  
  // Specify which fields will advance in time
  simulationState.SetStateFields("nodeScalar nodeVector0 nodeVector1");
  size_t numBufferNodes = pbsGrid.BufferSize();
  simulationState.Create(numBufferNodes,0);

  std::vector<size_t> partitionSizes(pbsPartInterval.Sizes());
  std::vector<size_t> partitionStarts(pbsPartInterval.Starts());
  std::vector<size_t> partitionBufferStarts(localPartitionInterval.Starts());
  
  testState = simulationState;


  // === Set up the initial state ===
  size_t iNode = 1024;
  double meshScalar = 12.0;
  double simScalar  = 15.0;
  std::vector<double> meshVector0(3,24.0);
  std::vector<double> meshVector1(3,1.0);
  std::vector<double> nodeScalar(numBufferNodes,0.0);
  std::vector<double> nodeVector0(3*numBufferNodes,0.0);
  std::vector<double> nodeVector1(3*numBufferNodes,0.0);
  std::vector<int>    nodeFlag(numBufferNodes,0);
  std::vector<double> simVector0(10,10.0);
  std::vector<double> simVector1(2,1.0);
  std::vector<char>   simBytes(10,'b');
  std::vector<int>    simInts(5,13);


  for(size_t iZ = 0;iZ < partitionSizes[2];iZ++){
    size_t nodeZ   = iZ + partitionStarts[2];
    size_t bufferZ = iZ + partitionBufferStarts[2];
    for(size_t iY = 0;iY < partitionSizes[1];iY++){
      size_t nodeY   = iY + partitionStarts[1];
      size_t bufferY = iY + partitionBufferStarts[1];
      for(size_t iX = 0;iX < partitionSizes[0];iX++){
        size_t nodeX   = iX + partitionStarts[0];
        size_t bufferX = iX + partitionBufferStarts[0];
        size_t bufferIndex = bufferZ*(bufferSizes[0]*bufferSizes[1]) +
          bufferY*bufferSizes[0] + bufferX;
        nodeScalar[bufferIndex]                    = 27.0;
        nodeVector0[bufferIndex]                   = 2.0;
        nodeVector0[bufferIndex+numPointsBuffer]   = 4.0;
        nodeVector0[bufferIndex+2*numPointsBuffer] = 8.0;
        nodeVector1[bufferIndex]                   = (myRank+1);
        nodeVector1[bufferIndex+numPointsBuffer]   = (myRank+1)*2.0;
        nodeVector1[bufferIndex+2*numPointsBuffer] = (myRank+1)*4.0;
        nodeFlag[bufferIndex]                      = myRank;
      }
    }
  }
  numComponents = 7;

  halo_t &testHalo(pbsGrid.Halo());
  int nodeDataMessageId = testHalo.CreateMessageBuffers(numComponents);

  testHalo.PackMessageBuffers(nodeDataMessageId,0,1,&nodeScalar[0]);
  testHalo.PackMessageBuffers(nodeDataMessageId,1,3,&nodeVector0[0]);
  testHalo.PackMessageBuffers(nodeDataMessageId,4,3,&nodeVector1[0]);

  testHalo.SendMessage(nodeDataMessageId,cartNeighbors,pbsGridComm);
  testHalo.ReceiveMessage(nodeDataMessageId,cartNeighbors,pbsGridComm);
  testHalo.UnPackMessageBuffers(nodeDataMessageId,0,1,&nodeScalar[0]);
  testHalo.UnPackMessageBuffers(nodeDataMessageId,1,3,&nodeVector0[0]);
  testHalo.UnPackMessageBuffers(nodeDataMessageId,4,3,&nodeVector1[0]);
  
  // testing data
  size_t tiNode = 0;
  double tmeshScalar = 0.0;
  double tsimScalar  = 0.0;
  std::vector<double> tmeshVector0(3,0.0);
  std::vector<double> tmeshVector1(3,0.0);
  std::vector<double> tnodeScalar(numBufferNodes,0.0);
  std::vector<double> tnodeVector0(3*numBufferNodes,0.0);
  std::vector<double> tnodeVector1(3*numBufferNodes,0.0);
  std::vector<int>    tnodeFlag(numBufferNodes,0);
  std::vector<double> tsimVector0(10,0);
  std::vector<double> tsimVector1(2,0);
  std::vector<char>   tsimBytes(10,0);
  std::vector<int>    tsimInts(5,0);

  // Link the domain's state storage with local storage
  simulationState.SetFieldBuffer("meshScalar",&meshScalar);
  simulationState.SetFieldBuffer("simScalar",&simScalar);
  simulationState.SetFieldBuffer("meshVector0",&meshVector0[0]);
  simulationState.SetFieldBuffer("meshVector1",&meshVector1[0]);
  simulationState.SetFieldBuffer("nodeScalar",&nodeScalar[0]);
  simulationState.SetFieldBuffer("nodeVector0",&nodeVector0[0]);
  simulationState.SetFieldBuffer("nodeVector1",&nodeVector1[0]);
  simulationState.SetFieldBuffer("nodeFlag",&nodeFlag[0]);
  simulationState.SetFieldBuffer("simVector0",&simVector0[0]);
  simulationState.SetFieldBuffer("simVector1",&simVector1[0]);
  simulationState.SetFieldBuffer("simInts",&simInts[0]);
  simulationState.SetFieldBuffer("simBytes",&simBytes[0]);

  testState.SetFieldBuffer("meshScalar",&tmeshScalar);
  testState.SetFieldBuffer("simScalar",&tsimScalar);
  testState.SetFieldBuffer("meshVector0",&tmeshVector0[0]);
  testState.SetFieldBuffer("meshVector1",&tmeshVector1[0]);
  testState.SetFieldBuffer("nodeScalar",&tnodeScalar[0]);
  testState.SetFieldBuffer("nodeVector0",&tnodeVector0[0]);
  testState.SetFieldBuffer("nodeVector1",&tnodeVector1[0]);
  testState.SetFieldBuffer("nodeFlag",&tnodeFlag[0]);
  testState.SetFieldBuffer("simVector0",&tsimVector0[0]);
  testState.SetFieldBuffer("simVector1",&tsimVector1[0]);
  testState.SetFieldBuffer("simInts",&tsimInts[0]);
  testState.SetFieldBuffer("simBytes",&tsimBytes[0]);

  std::ostringstream stateMetaDataStream;
  pcpp::field::metadataset &stateMetaData(simulationState.Meta());
  stateMetaData.WriteToStream(stateMetaDataStream);
  std::string stateMetaDataString(stateMetaDataStream.str());

  if(ix::sys::FILEEXISTS("statemetadata.h5"))
    ix::sys::Remove("statemetadata.h5");
  pcpp::io::hdf5::base metaFile("statemetadata.h5",pbsGridComm);
  if(!stateMetaDataString.empty())
    metaFile.CreateAttribute("StateMetaData",stateMetaDataString);
  metaFile.Close();

  state_t inputState;
  metaFile.Open("statemetadata.h5",pbsGridComm,true);
  std::string inputMetaDataString;
  metaFile.ReadAttribute("StateMetaData",inputMetaDataString);
  std::istringstream stateMetaDataInStream(inputMetaDataString);
  pcpp::field::metadataset &inputMetaData(inputState.Meta());
  inputMetaData.ReadFromStream(stateMetaDataInStream);
  metaFile.Close();

  bool metaDataIO = true;
  std::ostringstream inputMetaDataOutStream;
  inputMetaData.WriteToStream(inputMetaDataOutStream);
  std::string inputMetaString(inputMetaDataOutStream.str());
  if(inputMetaData != stateMetaData){
    std::cout << "MetaData: " << inputMetaString << std::endl;
    metaDataIO = false;
  }  
  metaDataIO = CheckResult(metaDataIO,testComm);
  parallelUnitResults.UpdateResult("HDF5:Interface:StateMetaDataIO",metaDataIO);

  if(ix::sys::FILEEXISTS("statedata.h5"))
    ix::sys::Remove("statedata.h5");
  pcpp::io::hdf5::base stateFile("statedata.h5",pbsGridComm);
  stateFile.CreateGroup("/StateData");
  stateFile.CreateGroup("/MeshData");
  pcpp::io::hdf5::WriteStateAttributes(simulationState,'s',"/StateData",stateFile);
  pcpp::io::hdf5::WriteStateAttributes(simulationState,'m',"/MeshData",stateFile);
  stateFile.Close();
  
  stateFile.Open("statedata.h5",pbsGridComm);
  pcpp::io::hdf5::ReadStateAttributes(testState,'s',"/StateData",stateFile);
  pcpp::io::hdf5::ReadStateAttributes(testState,'m',"/MeshData",stateFile);
  stateFile.Close();
  
  bool stateAttributeIO = true;
  if(tsimVector0 != simVector0)
    stateAttributeIO = false;
  if(tsimVector1 != simVector1)
    stateAttributeIO = false;
  if(tsimScalar != simScalar)
    stateAttributeIO = false;
  if(tsimBytes != simBytes)
    stateAttributeIO = false;
  if(tsimInts  != simInts)
    stateAttributeIO = false;
  if(tmeshScalar != meshScalar)
    stateAttributeIO = false;
  if(tmeshVector0 != meshVector0)
    stateAttributeIO = false;
  if(tmeshVector1 != meshVector1)
    stateAttributeIO = false;

  stateAttributeIO = CheckResult(stateAttributeIO,testComm);
  parallelUnitResults.UpdateResult("HDF5:Interface:StateAttributeIO",stateAttributeIO);
  
  
  if(ix::sys::FILEEXISTS("statenodedata.h5"))
    ix::sys::Remove("statenodedata.h5");
  pcpp::io::hdf5::base nodeFile("statenodedata.h5",pbsGridComm);
  nodeFile.CreateGroup("/MeshData");

  pcpp::io::hdf5::WriteStateNodeData(simulationState,gridSize,partitionSizes,partitionStarts,
                                     bufferSizes,partitionBufferStarts,"/MeshData",nodeFile);
  nodeFile.Close();
  
  
  nodeFile.Open("statenodedata.h5",pbsGridComm);
  pcpp::io::hdf5::ReadStateNodeData(testState,gridSize,partitionSizes,partitionStarts,
                                    bufferSizes,partitionBufferStarts,"/MeshData",nodeFile);
  nodeFile.Close();

  testHalo.PackMessageBuffers(nodeDataMessageId,0,1,&tnodeScalar[0]);
  testHalo.PackMessageBuffers(nodeDataMessageId,1,3,&tnodeVector0[0]);
  testHalo.PackMessageBuffers(nodeDataMessageId,4,3,&tnodeVector1[0]);

  testHalo.SendMessage(nodeDataMessageId,cartNeighbors,pbsGridComm);
  testHalo.ReceiveMessage(nodeDataMessageId,cartNeighbors,pbsGridComm);
  testHalo.UnPackMessageBuffers(nodeDataMessageId,0,1,&tnodeScalar[0]);
  testHalo.UnPackMessageBuffers(nodeDataMessageId,1,3,&tnodeVector0[0]);
  testHalo.UnPackMessageBuffers(nodeDataMessageId,4,3,&tnodeVector1[0]);

  bool stateNodeDataIO = true;
  if(tnodeScalar != nodeScalar)
    stateNodeDataIO = false;
  if(tnodeVector0 != nodeVector0)
    stateNodeDataIO = false;
  if(tnodeVector1 != nodeVector1)
    stateNodeDataIO = false;
  if(tnodeFlag != nodeFlag)
    stateNodeDataIO = false;
  

  stateNodeDataIO = CheckResult(stateNodeDataIO,testComm);
  parallelUnitResults.UpdateResult("HDF5:Interface:StateNodeDataIO",stateNodeDataIO);

  if(myRank == 0){
    if(ix::sys::FILEEXISTS("pbsgridtest.h5"))
      ix::sys::Remove("pbsgridtest.h5");
    if(ix::sys::FILEEXISTS("statenodedata.h5"))
      ix::sys::Remove("statenodedata.h5");
  }
  myGlobal.Finalize();
}

void TestPlasCom2Output(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{

  std::ostringstream messageStream;
  bool pc2OutputTest = true;

  typedef simulation::grid::parallel_blockstructured grid_t;
  typedef plascom2::operators::sbp::base             operator_t;
  typedef simulation::grid::halo                     halo_t;
  typedef simulation::state::base                    state_t;
  typedef pcpp::CommunicatorType                     comm_t;
  typedef pcpp::IndexIntervalType                    interval_t;
  typedef fixtures::ConfigurationType                config_t;

  typedef simulation::domain::base<grid_t,state_t>  domain_t;
  typedef euler::rhs<grid_t,state_t,operator_t>     rhs_t;
  
  config_t testConfig;
  config_t domainConfig;
  config_t geomConfig;

  global_t myGlobal;

  myGlobal.Init("TestPC2Output",testComm);
  myGlobal.SetVerbLevel(3);
  myGlobal.Profiling(true);

  int numProc = testComm.Size();
  int myRank  = testComm.Rank();


  std::ostringstream configMessageStream;
  if(pcpp::util::ProcessConfigurationFile("Config/Examples/misc/plascom2.test.config",testConfig,testComm,configMessageStream)){
    pc2OutputTest = false;
  }
  pc2OutputTest = CheckResult(pc2OutputTest,testComm);
  if(!pc2OutputTest){
    std::cout << "TestPlasCom2Output:ERROR!! ProcessConfigurationFile failed. Skipping tests." 
              << std::endl;
    parallelUnitResults.UpdateResult("HDF5:Output",pc2OutputTest);
    if(myRank == 0){
      std::cout << "Failed processing Config/Examples/misc/plascom2.test.config: " 
                << std::endl
                << configMessageStream.str() << std::endl << "Skipping remaining tests." << std::endl;
    }
    return;
  }
  //   for(int iRank = 0;iRank < numProc;iRank++){
  //     if(iRank == myRank){
  //       std::cout << " ======= Rank " << myRank << " =========" << std::endl
  //                 << configMessageStream.str() << std::endl;
  //     }
  //     testComm.Barrier();
  //   }
  
  std::vector<std::string> domainNames = testConfig.GetValueVector<std::string>("Domain:Names");
  std::vector<std::string> geometryNames = testConfig.GetValueVector<std::string>("Geometry:Names");

  int numDomains = domainNames.size();
  int numGeometries = geometryNames.size();

  for(int iDomain = 0;iDomain < numDomains;iDomain++){
    std::string &domainName(domainNames[iDomain]);
    domainConfig = pcpp::util::GetSubConfig(domainName,testConfig);
    if(myRank == 0){
      std::cout << "SubConfig(" << domainName << "):" << std::endl;
      std::cout << domainConfig << std::endl;
    }
  }
  int numDomainGrids = 1;

  std::string domainName(domainNames[0]);
  std::string domainGridName(domainConfig.GetValue<std::string>("GridNames"));

  for(int iGeom = 0;iGeom < numGeometries;iGeom++){
    std::string &geometryName(geometryNames[iGeom]);
    geomConfig = pcpp::util::GetSubConfig(geometryName,testConfig);
    if(myRank == 0){
      std::cout << "SubConfig(" << geometryName << "):" << std::endl;
      std::cout << geomConfig << std::endl;
    }
  }
  if(geometryNames.empty()){
    geometryNames.push_back("test-geometry");
  }
  std::string geometryName(geometryNames[0]);
  std::vector<std::string> gridNames = geomConfig.GetValueVector<std::string>("GridNames");
  if(gridNames.empty()){
    pc2OutputTest = false;
  }
  pc2OutputTest = CheckResult(pc2OutputTest,testComm);
  if(!pc2OutputTest){
    std::cout << "TestPlasCom2Output:ERROR!! Get gridnames fails. Skipping tests."
              << std::endl;
    parallelUnitResults.UpdateResult("HDF5:Output",pc2OutputTest);
    if(myRank == 0)
      std::cout << "Failed reading grid names, exiting" << std::endl;
    return;
  }
  
  std::string gridName(gridNames[0]);
  config_t gridConfig(pcpp::util::GetSubConfig(gridName,geomConfig));
  
  int numDim = gridConfig.GetValue<double>("Dimension");
  std::vector<double> physicalExtent(gridConfig.GetValueVector<double>("PhysicalExtent"));
  if(myRank == 0) {
    std::cout << "Grid Physical Extent:(";
    pcpp::io::DumpContents(std::cout,physicalExtent,",");
    std::cout << ")" << std::endl;
  }
  
  std::vector<size_t> gridSizes(gridConfig.GetValueVector<size_t>("NumPoints"));

  //  std::vector<double> dX(numDim,0);
  //   size_t numGlobalNodes = 1;
  //   size_t numGlobalCells = 1;
  //   for(int iDim = 0;iDim < numDim;iDim++){
  //     dX[iDim] = (physicalExtent[iDim*2+1] - physicalExtent[iDim*2])/static_cast<double>(gridSizes[iDim] - 1);
  //     numGlobalNodes *= gridSizes[iDim]; 
  //     numGlobalCells *= (gridSizes[iDim]-1);
  //   }
  
  grid_t testGrid;
  testGrid.SetGridSizes(gridSizes);
  testGrid.SetType(simulation::grid::UNIRECT);

  int operatorOrder = domainConfig.GetValue<int>("SpatialOrder");
  operator_t sbpOperator;
  plascom2::operators::sbp::Initialize(sbpOperator,operatorOrder);
  int boundaryDepth = (sbpOperator.numStencils-1)/2;
  std::vector<int> haloDepths(2*numDim,boundaryDepth);

  
  // Set up parallel Cartesian topology
  pcpp::ParallelTopologyInfoType cartInfo;
  cartInfo.numDimensions = numDim;
  cartInfo.isPeriodic = domainConfig.GetValueVector<int>("PeriodicDirs");
  testGrid.ParallelSetup(testComm,cartInfo,boundaryDepth,messageStream);

  comm_t &gridComm(testGrid.Communicator());
  
  // Grab the local coordinates and global dimension of the Cartesian topo
  std::vector<int> &cartCoords(gridComm.CartCoordinates());
  std::vector<int> &cartDims(gridComm.CartDimensions());

  // Grab the +/- neighbors in each Cartesian direction
  std::vector<int> cartNeighbors;
  gridComm.CartNeighbors(cartNeighbors);

  // Generate a report of the Cartesian setup and put it to stdout
  pcpp::report::CartesianSetup(messageStream,cartInfo);
  myGlobal.StdOut(messageStream.str(),2);
  pcpp::io::RenewStream(messageStream);

  interval_t globalInterval;
  globalInterval.InitSimple(gridSizes);

  interval_t &partInterval(testGrid.PartitionInterval());

  std::ostringstream gridStream;
  testGrid.GenerateCoordinates(gridStream);

  
  const std::vector<size_t> &bufferSizes(testGrid.BufferSizes());
  const pcpp::IndexIntervalType &partitionBufferInterval(testGrid.PartitionBufferInterval());  

  size_t numPointsBuffer = testGrid.BufferSize();

  // Set up the HALO data structure 
  halo_t &testHalo(testGrid.Halo());
 
  
  const std::vector<pcpp::IndexIntervalType> &remoteHaloBufferExtents(testHalo.RemoteHaloBufferExtents());
  const std::vector<pcpp::IndexIntervalType> &localHaloBufferExtents(testHalo.LocalHaloBufferExtents());

  pcpp::IndexIntervalType bufferInterval;
  bufferInterval.InitSimple(bufferSizes);

//   const std::vector<std::vector<size_t> > &recvIndices(testHalo.RecvIndices());

  // Create simple message buffers - assuming the grid coordinates will
  // be used to test 
  // int numComponents = numDim+2;
  // int xyzMessageId = testHalo.CreateMessageBuffers(numComponents);

  //  testGrid.AllocateCoordinateData();
  std::vector<double> &gridCoordinates(testGrid.CoordinateData());
  std::vector<double> &dX(testGrid.DX());
  
  
  state_t    simState;
  state_t    paramState;

  if(testfixtures::euler::SetupEulerState(testGrid,simState,paramState,0)){
    pc2OutputTest = false;
  }
  pc2OutputTest = CheckResult(pc2OutputTest,testComm);
  if(!pc2OutputTest){
    std::cout << "TestPlasCom2Output:ERROR!! SetupEulerState failed. Skipping tests." 
              << std::endl;
    parallelUnitResults.UpdateResult("HDF5:Output",pc2OutputTest);
    if(myRank == 0)
      std::cout << "Failed Euler state setup, exiting" << std::endl;
    return;
  }

  double t = 0.0;
  double gamma = 1.4;
  double dt = .001;
  double cfl = 1.0;
  double Re  = 0.0;
  int flags = 1;

  // Add some fields that are not usual
  paramState.AddField("dX",'s',3,8,"space");
  paramState.Create(numPointsBuffer,0);

  paramState.SetFieldBuffer("gamma",&gamma);
  paramState.SetFieldBuffer("inputDT",&dt);
  paramState.SetFieldBuffer("inputCFL",&cfl);
  paramState.SetFieldBuffer("Flags",&flags);  
  paramState.SetFieldBuffer("refRe",&Re);
  paramState.SetFieldBuffer("dX",&dX[0]);

  std::vector<double> rho(numPointsBuffer,1.0);
  std::vector<double> rhoV(3*numPointsBuffer,0.0);
  std::vector<double> rhoE(3*numPointsBuffer,1.0/(gamma-1.0));

  std::vector<double> dRho(numPointsBuffer,1.0);
  std::vector<double> dRhoV(3*numPointsBuffer,0.0);
  std::vector<double> dRhoE(3*numPointsBuffer,1.0/(gamma-1.0));

  std::vector<double> T(numPointsBuffer,0.0);
  std::vector<double> p(numPointsBuffer,0.0);
  std::vector<double> V(3*numPointsBuffer,0.0);

  std::vector<double> rhom1(numPointsBuffer,0.0);


  simState.SetFieldBuffer("simTime",&t);
  simState.SetFieldBuffer("rho",&rho[0]);
  simState.SetFieldBuffer("rhoV",&rhoV[0]);
  simState.SetFieldBuffer("rhoE",&rhoE[0]);

  simState.SetFieldBuffer("pressure",&p[0]);
  simState.SetFieldBuffer("temperature",&T[0]);
  simState.SetFieldBuffer("velocity",&V[0]);
  simState.SetFieldBuffer("rhom1",&rhom1[0]);

  state_t rhsState(simState);

  rhsState.SetFieldBuffer("rho",&dRho[0]);
  rhsState.SetFieldBuffer("rhoV",&dRhoV[0]);
  rhsState.SetFieldBuffer("rhoE",&dRhoE[0]);

  // set up the domain - don't care about this in testing
  //   domain_t eulerDomain;
  //   eulerDomain.numGrids = 1;
  //   eulerDomain.Grids().resize(1,&simGrid);
  //   eulerDomain.States().resize(1,&simState);

  // -- Set up the RHS --
  rhs_t eulerRHS;
  
  if(eulerRHS.Initialize(testGrid,simState,paramState,sbpOperator)){
    pc2OutputTest = false;
  }

  pc2OutputTest = CheckResult(pc2OutputTest,testComm);
  if(!pc2OutputTest){
    std::cout << "TestPlasCom2Output:ERROR!! Failed EulerRHS::Initialize failed. Skipping tests."
              << std::endl;
    parallelUnitResults.UpdateResult("HDF5:Output",pc2OutputTest);
    if(myRank == 0)
      std::cout << "Failed Euler initialization, exiting" << std::endl;
    return;
  }  
  if(eulerRHS.SetupRHS(testGrid,simState,rhsState)){
    pc2OutputTest = false;
  }
  pc2OutputTest = CheckResult(pc2OutputTest,testComm);
  if(!pc2OutputTest){
    std::cout << "TestPlasCom2Output:ERROR!! Failed EulerRHS::SetupRHS failed. Skipping tests."
              << std::endl;
    parallelUnitResults.UpdateResult("HDF5:Output",pc2OutputTest);
    if(myRank == 0)
      std::cout << "Failed Euler RHS setup, exiting" << std::endl;
    return;
  }  
  std::cout << "TestPlasCom2Output: Getting Euler RHS" << std::endl;
  if(eulerRHS.RHS(0)){
    pc2OutputTest = false;
  }
  std::cout << "TestPlasCom2Output: Euler RHS returned." << std::endl;
  pc2OutputTest = CheckResult(pc2OutputTest,testComm);
  if(!pc2OutputTest){
    std::cout << "TestPlasCom2Output:ERROR!! Failed EulerRHS::RHS failed. Skipping tests."
              << std::endl;
    parallelUnitResults.UpdateResult("HDF5:Output",pc2OutputTest);
    if(myRank == 0)
      std::cout << "Failed Euler RHS setup 2, exiting" << std::endl;
    return;
  }  
  
  // Begin PlasCom2 output
  if(myRank == 0){
    if(ix::sys::FILEEXISTS("plascom2_output.h5"))
      ix::sys::Remove("plascom2_output.h5");
  }
  testComm.Barrier();

  pcpp::io::hdf5::base pc2File("plascom2_output.h5",gridComm);
  pc2File.CreateGroup("/PlasCom2");

  pc2File.CreateGroup("/PlasCom2/Configuration");
  std::ostringstream configStream;
  configStream << testConfig;
  std::string configString(configStream.str());

  if(!configString.empty()){
    //     if(testComm.Rank() == 0){
    //       std::cout << "writing config string: ((((" << std::endl
    //                 << configString << std::endl << "))))" << std::endl;
    //       std::cout << "COnfig size: " << configString.size() << std::endl;
    //     }
    pc2File.CreateAttribute("/PlasCom2/Configuration/Parameters",configString);
  }
  
  pc2File.CreateGroup("/PlasCom2/Geometry");
  std::vector<size_t> attSize(1,1);
  pc2File.CreateAttribute("/PlasCom2/Geometry/numGeometries",attSize,&numGeometries);
  pc2File.CreateAttribute("/PlasCom2/Geometry/geometryNames",geometryName);
  std::string geomGroupName(std::string("/PlasCom2/Geometry/")+geometryName);
  pc2File.CreateGroup(geomGroupName);
  int numGrids = 1;
  pc2File.CreateAttribute(geomGroupName+"/numGrids",attSize,&numGrids);
  pc2File.CreateAttribute(geomGroupName+"/gridNames",gridName);
  std::vector<double> expectedCoordinates(testGrid.CoordinateData());
  pcpp::io::hdf5::WriteGrid(testGrid,gridName,geomGroupName,pc2File);

  pc2File.CreateGroup("/PlasCom2/Simulation");
  pc2File.CreateAttribute("/PlasCom2/Simulation/numDomains",attSize,&numDomains);
  pc2File.CreateAttribute("/PlasCom2/Simulation/domainNames",domainName);
  pcpp::io::hdf5::WriteStateAttributes(simState,'s',"/PlasCom2/Simulation",pc2File);
  pcpp::io::hdf5::WriteStateAttributes(paramState,'s',"/PlasCom2/Simulation",pc2File);
  std::string domainGroupName(std::string("/PlasCom2/Simulation/")+domainName);
  pc2File.CreateGroup(domainGroupName);
  pc2File.CreateGroup(domainGroupName+"/Configuration");
  std::ostringstream domainConfigStream;
  domainConfigStream << domainConfig;
  std::string domainConfigString(domainConfigStream.str());
  if(!domainConfigString.empty())
    pc2File.CreateAttribute(domainGroupName+"/Configuration/Parameters",domainConfigString);
  pc2File.CreateAttribute(domainGroupName+"/numGrids",attSize,&numDomainGrids);
  pc2File.CreateAttribute(domainGroupName+"/gridNames",domainGridName);
  pcpp::io::hdf5::WriteStateAttributes(simState,'d',domainGroupName,pc2File);
  pcpp::io::hdf5::WriteStateAttributes(paramState,'d',domainGroupName,pc2File);
  std::string gridGroupName(domainGroupName+"/"+domainGridName);
  pc2File.CreateGroup(gridGroupName);
  std::vector<size_t> partitionSizes(partInterval.Sizes());
  std::vector<size_t> partitionStarts(partInterval.Starts());
  std::vector<size_t> partitionBufferStarts(partitionBufferInterval.Starts());
  pcpp::io::hdf5::WriteStateNodeData(simState,gridSizes,partitionSizes,partitionStarts,
                                     bufferSizes,partitionBufferStarts,gridGroupName,pc2File);  
  pc2File.Close();
  

  
  pc2OutputTest = CheckResult(pc2OutputTest,testComm);
  if(!pc2OutputTest){
    if(myRank == 0)
      std::cout << "Failed PlasCom2 output, exiting" << std::endl;
  }
  parallelUnitResults.UpdateResult("HDF5:Output",pc2OutputTest);

  if(myRank == 0){
    if(ix::sys::FILEEXISTS("plascom2_output.h5"))
      ix::sys::Remove("plascom2_output.h5");
  }

  return;
}

void TestPlasCom2Input(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{

  typedef simulation::grid::parallel_blockstructured grid_t;
  typedef plascom2::operators::sbp::base             operator_t;
  typedef simulation::grid::halo                     halo_t;
  typedef simulation::state::base                    state_t;
  typedef pcpp::CommunicatorType                     comm_t;
  typedef pcpp::IndexIntervalType                    interval_t;
  typedef fixtures::ConfigurationType                config_t;
  typedef simulation::domain::base<grid_t,state_t>   domain_t;
  typedef euler::rhs<grid_t,state_t,operator_t>      rhs_t;

  std::string functionName("TestPlasCom2Input");

  global_t myGlobal;
  myGlobal.Init(functionName,testComm);
  myGlobal.SetVerbLevel(3);
  myGlobal.Profiling(true);

  int myRank = testComm.Rank();

  myGlobal.StdOut(std::string("TestFunction: ")+functionName+std::string("\n"));
  
  std::string originalDirectory(ix::sys::CWD());

  ix::sys::ChDir("InviscidTests/AcousticPulseFarfield2D/PlasCom2/");

  testComm.Barrier();

  bool processConfigurationRestart = true;

  std::ostringstream messageStream;
  std::string configFileName("curvilinear.config");
  pcpp::ConfigType referenceConfig;

  if(ProcessConfigurationFile(configFileName,referenceConfig,testComm,messageStream)){
    myGlobal.StdOut(std::string("Failed to process reference configuration from curvilinear.config.\nMessages:")+
                    messageStream.str()+std::string("\n"));
    processConfigurationRestart = false;
  } 

  if(myRank == 0){

    if(ix::sys::FILEEXISTS("tmptest.h5"))
      ix::sys::Remove("tmptest.h5");

    pcpp::io::hdf5::base testFile("tmptest.h5");
    testFile.CreateGroup("/PlasCom2");
    testFile.CreateGroup("/PlasCom2/Geometry");
    testFile.CreateGroup("/PlasCom2/Domains");
    testFile.CreateGroup("/PlasCom2/Geometry/Grids");
    testFile.CreateGroup("/PlasCom2/Configuration");
    std::ostringstream configStream;
    configStream << referenceConfig;
    std::string configString(configStream.str());
    if(!configString.empty())
      testFile.CreateAttribute("/PlasCom2/Configuration/Parameters",configString);
    testFile.Close();
  }

  testComm.Barrier();

  pcpp::ConfigType pc2Config;
  std::string testingFileName("tmptest.h5");
  if(!ix::sys::FILEEXISTS(testingFileName)){
    myGlobal.StdOut(std::string("Could not find testing file (")+testingFileName+
                    std::string(").\n"));
    processConfigurationRestart = false;
  }
  if(pcpp::util::ProcessConfigurationRestart(testingFileName,pc2Config,testComm,messageStream)){
    myGlobal.StdOut(std::string("Failed to process restart file (")+testingFileName+
                    std::string(").\nMessages:\n")+messageStream.str()+std::string("\n"));
    processConfigurationRestart = false;
  }
  pcpp::io::RenewStream(messageStream);

  if(processConfigurationRestart){
    if(pc2Config != referenceConfig){
      myGlobal.StdOut("Configuration did not match reference configuration. Test failed.\n");
      processConfigurationRestart = false;
      std::ostringstream Ostr;
      Ostr << pc2Config << std::endl;
      myGlobal.StdOut(Ostr.str());
    }
  }

  parallelUnitResults.UpdateResult("HDF5:Restart:Configuration",processConfigurationRestart);

  testComm.Barrier();

  bool fileInfoStreamTest = true;
  pcpp::io::simfileinfo restartFileInfo;
  std::string restartFileName("accepted/PlasCom2_000001000.h5");
  int statusCode = 0;
  
  if(myRank == 0){
    std::ostringstream errStream;
    if(pcpp::io::hdf5::FileInfo(restartFileName,restartFileInfo,errStream)){
      messageStream << "ERROR: FileInfo failed for " << restartFileName 
                    << ".\nHDF5 messages:" 
                    << errStream.str() << std::endl;
      myGlobal.StdOut(messageStream.str());
      pcpp::io::RenewStream(messageStream);
      testComm.SetErr(1);
    } else {
      std::ostringstream infoOut;
      infoOut << "Original info for " << restartFileName << ": " << std::endl;
      pcpp::report::SimFileInfo(infoOut,restartFileInfo);
      infoOut << "RAW stream output for info:" << std::endl;
      std::ostringstream infoStream;
      infoStream << restartFileInfo;
      infoOut << infoStream.str() << std::endl
              << "------------------" << std::endl;
      myGlobal.StdOut(infoOut.str());
      std::istringstream infoInStream(infoStream.str());
      pcpp::io::simfileinfo testInfo;
      infoInStream >> testInfo;
      std::ostringstream testStream;
      testStream << testInfo;
      if(!(testInfo == restartFileInfo)){
        testComm.SetErr(1);
        myGlobal.StdOut("File info object failed stream tests.\n");
      } else {
        myGlobal.StdOut("RAW stream output for test:\n");
        myGlobal.StdOut(testStream.str()+std::string("\n"));
        myGlobal.StdOut("=========\n");
      }
    } 
  }
  
  if(testComm.Check()){
    statusCode = 1;
    fileInfoStreamTest = false;
    testComm.SetErr(0);
  }

  parallelUnitResults.UpdateResult("HDF5:Restart:FileInfoStream",fileInfoStreamTest);

  bool fileInfoBroadcastTest = true;
  pcpp::io::simfileinfo myFileInfo;

  if(fileInfoStreamTest){
    if(!ix::sys::FILEEXISTS(restartFileName)){
      myGlobal.ErrOut(std::string("Could not find required restart file: ")+restartFileName+
                      std::string("\n"));
      
      fileInfoBroadcastTest = false;
    } else {
      std::ostringstream errStream;
      testComm.StreamBroadCast(restartFileInfo);
      
      int errCode = pcpp::io::hdf5::FileInfo(restartFileName,myFileInfo,testComm,errStream);

      if(errCode > 0){
        std::cout << "FileInfo failed on rank " << myRank << std::endl;
      } else {
        if(!(restartFileInfo == myFileInfo)){
          std::cout << "My file info and broadcasted one do not match." << std::endl;
          testComm.SetErr(1);
        } else {
          std::cout << "info broadcast succeeded on rank " << myRank << std::endl;
        }
        if(testComm.Check())
          fileInfoBroadcastTest = false;
      }
    }

    if(!fileInfoBroadcastTest){
      std::string testResultString;
      if(myRank == 1){
        std::ostringstream infoOut;
        infoOut << "Original: " << std::endl;
        pcpp::report::SimFileInfo(infoOut,myFileInfo);
        infoOut << std::endl << "Broadcast version: " << std::endl;
        pcpp::report::SimFileInfo(infoOut,restartFileInfo);
        infoOut << std::endl;
        infoOut << "Original RAW stream output:" << std::endl
                << "--------" << std::endl
                << myFileInfo << std::endl
                << "+++++++"  << std::endl
                << "BC RAW stream output: " << std::endl
                << "------"   << std::endl
                << restartFileInfo << std::endl
                << "=======" << std::endl;
        testResultString = infoOut.str();    
      }
      testComm.BroadCast(testResultString,1);
      myGlobal.StdOut(testResultString);
    }
  } else {
    myGlobal.StdOut("File info broadcast test is skipped due to earlier errors.\n");
    fileInfoBroadcastTest = false;
  }

  parallelUnitResults.UpdateResult("HDF5:Restart:FileInfoBroadCast",fileInfoBroadcastTest);

  testComm.Barrier();

  ix::sys::ChDir(originalDirectory);

  bool readLegacyBCDat = true;
  const std::string bcFileName("PlasComCM/inputs/legacy.bc.dat");
  std::vector<plascom2::util::legacybc> bcData;
  if(plascom2::util::ReadLegacyBCDat(bcFileName,bcData,testComm,messageStream)){
    messageStream << "Read of legacy bc data failed." << std::endl;
    myGlobal.StdOut(messageStream.str());
    pcpp::io::RenewStream(messageStream);
    readLegacyBCDat = false;
  }
  
  std::ostringstream expectedBCDat;

  expectedBCDat 
    << "1    21     1     1     1     1    -1     1    -1" << std::endl
    << "1    22    -1    -1    -1     1    -1     1    -1" << std::endl
    << "1    24     2     1    -1     1     1     1    -1" << std::endl
    << "1    24    -2     1    -1    -1    -1     1    -1" << std::endl
    << "1    24     3     1    -1     1    -1     1     1" << std::endl
    << "1    24    -3     1    -1     1    -1    -1    -1" << std::endl
    << "2    21     1     1     1     1    -1     1    -1" << std::endl
    << "2    22    -1    -1    -1     1    -1     1    -1" << std::endl
    << "2    24     2     1    -1     1     1     1    -1" << std::endl
    << "2    24    -2     1    -1    -1    -1     1    -1" << std::endl
    << "2    24     3     1    -1     1    -1     1     1" << std::endl
    << "2    24    -3     1    -1     1    -1    -1    -1" << std::endl
    << "2    99     3     1    -1     1    -1     1    20" << std::endl
    << "2    99    -3     1    -1     1    -1    -20   -1" << std::endl;
  
  std::istringstream expectedBCStream(expectedBCDat.str());
  plascom2::util::legacybc myBcDat;
  std::vector<plascom2::util::legacybc> expectedBCData;
  while(expectedBCStream >> myBcDat)
    expectedBCData.push_back(myBcDat);
  if(bcData != expectedBCData)
    readLegacyBCDat = false;
  
  parallelUnitResults.UpdateResult("HDF5:Restart:ReadLegacyBCDat",readLegacyBCDat);
  
  ix::sys::ChDir("PlasComCM/examples/FlowPastCylinder");

  bool processLegacyConfigurationFile = true;
  std::string legacyConfigFileName("flowpastcylinder.inp");
  pcpp::ConfigType legacyConfig;
  if(pcpp::util::ProcessConfigurationFile(legacyConfigFileName,legacyConfig,testComm,messageStream)){
    myGlobal.StdOut(std::string("Failed to process legacy configuration from ")+
                    legacyConfigFileName+std::string(".\nMessages:")+
                    messageStream.str()+std::string("\n"));
    pcpp::io::RenewStream(messageStream);
    processLegacyConfigurationFile = false;
  }
  myGlobal.StdOut(messageStream.str());
  pcpp::io::RenewStream(messageStream);
  if(plascom2::util::ConvertLegacyConfiguration(legacyConfig,testComm,messageStream)){
    myGlobal.StdOut(std::string("Failed to convert legacy configuration from ")+
                    legacyConfigFileName+std::string(".\nMessages:")+
                    messageStream.str()+std::string("\n"));
    pcpp::io::RenewStream(messageStream);
    processLegacyConfigurationFile = false;
  }    
  messageStream << legacyConfig << std::endl;
  myGlobal.StdOut(messageStream.str());
  pcpp::io::RenewStream(messageStream);
 
  parallelUnitResults.UpdateResult("HDF5:Restart:ProcessLegacyConfig",processLegacyConfigurationFile);
  
  ix::sys::ChDir(originalDirectory);

  bool fetchGridData = true;
  restartFileName.assign("HDF5Examples/example_pc2_singlegrid.h5");
  pcpp::ConfigType exampleConfig;
  if(pcpp::util::ProcessConfigurationRestart(restartFileName,exampleConfig,testComm,messageStream)){
    myGlobal.ErrOut(messageStream.str()+std::string("\n"));
    myGlobal.ErrOut(std::string("Processing configuration from restart file (")+
                    restartFileName+std::string(") failed.\n"));
    testComm.SetErr(1);
  }
  if(testComm.Check())
    return;

  {
    bool parallelFileInfo = true;
    pcpp::io::simfileinfo singleGridFileInfo;
    if(pcpp::io::hdf5::FileInfo(restartFileName,singleGridFileInfo,testComm,messageStream))
      return;
    
    if(singleGridFileInfo.numGrids != 1)
      parallelFileInfo = false;
    const std::string &gridName(singleGridFileInfo.gridNames[0]);
    if(gridName != "box:grid3")
      parallelFileInfo = false;
    std::vector<size_t> expectedGridSizes(2,101);
    if(singleGridFileInfo.gridSizes[0] != expectedGridSizes)
      parallelFileInfo = false;
    
    parallelFileInfo = CheckResult(parallelFileInfo,testComm);
    
    parallelUnitResults.UpdateResult("HDF5:Restart:ParallelFileInfo",parallelFileInfo);

    grid_t testGrid;
    int haloDepth = 2;
    std::vector<size_t> gridSizes(exampleConfig.GetValueVector<size_t>("box:grid3:NumPoints"));
    int numDim = gridSizes.size();
    
    // Make a non-periodic cartesian grid & topology
    pcpp::ParallelTopologyInfoType cartInfo;
    cartInfo.numDimensions = 2;
    cartInfo.isPeriodic.resize(numDim,0);
    
    
    // Set grid sizes
    size_t numGlobalNodes = 1;
    size_t numGlobalCells = 1;
    for(int iDim = 0;iDim < numDim;iDim++){
      numGlobalNodes *= gridSizes[iDim];
      numGlobalCells *= (gridSizes[iDim]-1);
    }
    testGrid.SetGridSizes(gridSizes);
    
    messageStream << "Setting up grid with sizes: (";
    pcpp::io::DumpContents(messageStream,gridSizes,",");
    messageStream << ")" << std::endl;
    myGlobal.StdOut(messageStream.str());
    pcpp::io::RenewStream(messageStream);
    
    
    // This call does everything to set up the grid and required halos
    messageStream << "====== Grid Parallel Setup ========= " << std::endl;
    if(testGrid.ParallelSetup(testComm,cartInfo,haloDepth,messageStream)){
      myGlobal.ErrOut("Grid::ParallelSetup failed with messages:\n");
      myGlobal.ErrOut(messageStream.str());
      pcpp::io::RenewStream(messageStream);
      //    sectionResult = false;
    }
    messageStream << "==================================== " << std::endl;
    
    // Generate a report of the Cartesian grid and topo setup
    messageStream << "------- Grid & Topo Report --------- " << std::endl;
    pcpp::report::CartesianSetup(messageStream,cartInfo);    
    messageStream << simulation::report::Grid(testGrid)      << std::endl
                  << "------------------------------------ " << std::endl;
    myGlobal.StdOut(messageStream.str());
    pcpp::io::RenewStream(messageStream);
    
    
    if(pcpp::io::hdf5::FetchGridData(testGrid,gridName,restartFileName,messageStream)){
      myGlobal.ErrOut(std::string("Fetch grid data failed:\n")+messageStream.str()+std::string("\n"));
      testComm.SetErr(1);
    } else {
      myGlobal.StdOut(std::string("Read grid data for (")+gridName+
                      std::string(") from ")+restartFileName+std::string(".\n"),2);
      myGlobal.StdOut(messageStream.str(),2);
      pcpp::io::RenewStream(messageStream);
    }
    if(testComm.Check())
      fetchGridData = false;
    
    parallelUnitResults.UpdateResult("HDF5:Restart:FetchSingleGrid:PC2",fetchGridData);
  }
  
  {
    
    bool fetchMultipleGridData = true;
    restartFileName.assign("HDF5Examples/example_pc2_multiplegrids.h5");
    pcpp::ConfigType exampleConfig;
    if(pcpp::util::ProcessConfigurationRestart(restartFileName,exampleConfig,
                                               testComm,messageStream)){
      myGlobal.ErrOut(messageStream.str()+std::string("\n"));
      myGlobal.ErrOut(std::string("Processing configuration from restart file (")+
                      restartFileName+std::string(") failed.\n"));
      testComm.SetErr(1);
    }
    if(testComm.Check())
      return;
    
    bool parallelFileInfo = true;
    pcpp::io::simfileinfo multipleGridFileInfo;
    if(pcpp::io::hdf5::FileInfo(restartFileName,multipleGridFileInfo,testComm,messageStream))
      return;
    
    int numGrids = multipleGridFileInfo.numGrids;
    if(numGrids != 2) {
      parallelFileInfo = false;
    } else {
      const std::vector<std::string> &gridNames(multipleGridFileInfo.gridNames);
      if(gridNames[0] != "box:grid3" ||
         gridNames[1] != "box:grid4")
        parallelFileInfo = false;
    }

    const std::vector<std::string> &gridNames(multipleGridFileInfo.gridNames);
    std::vector<size_t> expectedGridSizes(2,101);

    for(int iGrid = 0;iGrid < numGrids;iGrid++)
      if(multipleGridFileInfo.gridSizes[iGrid] != expectedGridSizes)
        parallelFileInfo = false;
    
    parallelFileInfo = CheckResult(parallelFileInfo,testComm);
    
    parallelUnitResults.UpdateResult("HDF5:Restart:MultipleGridFileInfo",parallelFileInfo);
    
    for(int iGrid = 0;iGrid < numGrids;iGrid++){
      grid_t testGrid;
      int haloDepth = 2;
      const std::string &gridName(gridNames[iGrid]);
      std::ostringstream Ostr;
      Ostr << gridName << ":NumPoints";
      std::vector<size_t> gridSizes(exampleConfig.GetValueVector<size_t>(Ostr.str()));
      int numDim = gridSizes.size();
      
      // Make a non-periodic cartesian grid & topology
      pcpp::ParallelTopologyInfoType cartInfo;
      cartInfo.numDimensions = 2;
      cartInfo.isPeriodic.resize(numDim,0);
    
    
      // Set grid sizes
      size_t numGlobalNodes = 1;
      size_t numGlobalCells = 1;
      for(int iDim = 0;iDim < numDim;iDim++){
        numGlobalNodes *= gridSizes[iDim];
        numGlobalCells *= (gridSizes[iDim]-1);
      }
      testGrid.SetGridSizes(gridSizes);
    
      messageStream << "Setting up grid with sizes: (";
      pcpp::io::DumpContents(messageStream,gridSizes,",");
      messageStream << ")" << std::endl;
      myGlobal.StdOut(messageStream.str());
      pcpp::io::RenewStream(messageStream);
      
      
      // This call does everything to set up the grid and required halos
      messageStream << "====== Grid Parallel Setup ========= " << std::endl;
      if(testGrid.ParallelSetup(testComm,cartInfo,haloDepth,messageStream)){
        myGlobal.ErrOut("Grid::ParallelSetup failed with messages:\n");
        myGlobal.ErrOut(messageStream.str());
        pcpp::io::RenewStream(messageStream);
        //    sectionResult = false;
      }
      messageStream << "==================================== " << std::endl;
      
      // Generate a report of the Cartesian grid and topo setup
      messageStream << "------- Grid & Topo Report --------- " << std::endl;
      pcpp::report::CartesianSetup(messageStream,cartInfo);    
      messageStream << simulation::report::Grid(testGrid)      << std::endl
                    << "------------------------------------ " << std::endl;
      myGlobal.StdOut(messageStream.str());
      pcpp::io::RenewStream(messageStream);
      
      
      if(pcpp::io::hdf5::FetchGridData(testGrid,gridName,restartFileName,messageStream)){
        myGlobal.ErrOut(std::string("Fetch grid data failed:\n")+messageStream.str()+std::string("\n"));
        testComm.SetErr(1);
      } else {
        myGlobal.StdOut(std::string("Read grid data for (")+gridName+
                        std::string(") from ")+restartFileName+std::string(".\n"),2);
        myGlobal.StdOut(messageStream.str(),2);
        pcpp::io::RenewStream(messageStream);
      }
      if(testComm.Check())
        fetchMultipleGridData = false;
    }
    parallelUnitResults.UpdateResult("HDF5:Restart:FetchMultipleGrids:PC2",fetchMultipleGridData); 

  }
    ix::sys::ChDir(originalDirectory);
  return;
};