Domain.H

Go to the documentation of this file.

#ifndef __DOMAIN_H__

#include "PCPPTypes.H"
#include "PCPPUtil.H"
#include "Boundary.H"
#include "Stencil.H"
#include "ConfigUtil.H"
#include "NavierStokesBC.H"
#include "RHS.H"
#include "Geometry.H"

using namespace pcpp::util;
using namespace fixtures;

namespace simulation {

namespace domain {

#define CONSTANT_DT  0
#define CONSTANT_CFL 1
  

  struct info {
    info() : numDomains(0) {};
    int numDomains;
    std::vector<std::string> domainNames;
    int numLocalDomains;
    std::vector<int> globalDomainIndex;
  };

  struct partitioninfo {
    std::vector<int>    numProcsGrid;
    std::vector<size_t> numPointsProc;
    std::vector<int>    gridColor;
    std::vector<int>                        gridIndexOfLocalGrids;
  };

  template<typename GridT,typename StateT>
  class base {

    template<typename DomainT>
    friend class DomainInitializer;
    
  public:

    typedef GridT                            GridType;
    typedef StateT                           StateType;
    typedef typename GridType::HaloType      HaloType;

    // some fixed types, smells bad
    typedef boundary::base                      BoundaryType;
    typedef boundary::bc::base                  BCType;
    typedef typename StateType::MetaDataSetType DataDictionaryType;
    typedef plascom2::operators::sbp::base      OperatorType;

    //    typedef euler::rhs<GridType,StateType,OperatorType> RHSType;
    typedef simulation::rhs::domain_base<GridType,StateType,OperatorType> RHSType;

    //     typedef BoundaryT  BoundaryType;
    //     typedef RHST       RHSType;

    typedef std::vector<GridType*>     gridvector;
    typedef std::vector<StateType*>    statevector;
    typedef std::vector<RHSType*>      rhsvector;

    //    typedef simulation::state::halo    HaloType;

  public:

    base() : numGrids(0), rhsPtr(&dummyRHS), domainTime(0), stepNumber(0), domainRank(0),
             geometryInfoPtr(NULL)
    { };

    StateType   &State(int iGrid)  { return (*gridStates[iGrid]);  };
    StateType   &Target(int iGrid) { return (*gridTargets[iGrid]); };
    StateType   &Param(int iGrid)  { return (*gridParams[iGrid]);  };
    GridType    &Grid(int iGrid)   { return (*domainGrids[iGrid]); };
    RHSType     &RHS(int iGrid)    { return (*gridRHSs[iGrid]);    };

    gridvector  &Grids()          { return (domainGrids);         };
    statevector &States()         { return (gridStates);          };
    statevector &Targets()        { return (gridTargets);         };
    statevector &Params()         { return (gridParams);          };
    rhsvector   &RHSs()           { return (gridRHSs);            };

    std::vector<BoundaryType>  &GridBoundaries(int iGrid) { return (domainBoundaries[iGrid]); };
    std::vector<BoundaryType>  &DomainBoundary(int iGrid) { return (domainBoundaries[iGrid]); };
    std::vector<std::vector<BoundaryType> > &DomainBoundaries() { return (domainBoundaries); };
    
    const std::vector<std::string> &BCNames() const {return(bcNames); };
    const std::vector<std::string> &BoundaryNames() const {return(boundaryNames); };

    //    const std::vector<std::string> &BoundaryNames() const {return(boundaryNames); };

    void SetGrid(int gridIndex,GridType &inGrid)
    {
      assert(gridIndex >= 0 && gridIndex < domainGrids.size());
      domainGrids[gridIndex] = &inGrid;
    };

    void SetGridState(int gridIndex,StateType &inState)
    {
      assert(gridIndex >= 0 && gridIndex < domainGrids.size());
      gridStates[gridIndex] = &inState;
    };

    void SetGridTarget(int gridIndex,StateType &inState)
    {
      assert(gridIndex >= 0 && gridIndex < domainGrids.size());
      gridTargets[gridIndex] = &inState;
    };

    void SetGridParams(int gridIndex,StateType &paramState)
    {
      assert(gridIndex >= 0 && gridIndex < domainGrids.size());
      gridParams[gridIndex] = &paramState;
    };
    
    std::string Name() const { return(domainName); };

    void SetGlobal(pcpp::ParallelGlobalType &inGlobal)
    {
      globalPtr = &inGlobal;
    };

    int SetNumGrids(int numberOfGrids)
    {
      if(numberOfGrids <= 0)
        return(1);
      numGrids = numberOfGrids;
      ownGrid.resize(numGrids,false);
      ownState.resize(numGrids,false);
      ownRHS.resize(numGrids,false);
      domainGrids.resize(numGrids,NULL);
      gridStates.resize(numGrids,NULL);
      gridTargets.resize(numGrids,NULL);
      gridParams.resize(numGrids,NULL);
      gridRHSs.resize(numGrids,NULL);
      domainBoundaries.resize(numGrids);
      gridCommunicators.resize(numGrids,domainCommunicator);
      return(0);
    };

    int ConfigureDomain(fixtures::ConfigurationType &inConfig,std::ostream &inMessageStream)
    {

      SetMessageStream(inMessageStream);
      std::ostream &messageStream(*messageStreamPtr);

      domainConfig = inConfig;
      domainName = domainConfig.GetValue("ConfigName");

      if(domainName.empty()){
        messageStream << "Name of domain unspecified." << std::endl;
        return(1);
      }

      spatialOrder = domainConfig.GetValue<int>(ConfigKey(domainName,"SpatialOrder"));
      int interiorOrder  = (spatialOrder - 1)*2;
      plascom2::operators::sbp::Initialize(domainOperator,interiorOrder);
      
      if(domainConfig.IsSet(ConfigKey(domainName,"ScalarNames")))
        scalarNames = domainConfig.GetValueVector<std::string>(ConfigKey(domainName,"ScalarNames"));
      
      int numGridsConfigured = 0;
      gridNames = domainConfig.GetValueVector<std::string>(ConfigKey(domainName,"GridNames"));
      if(!domainConfig.IsSet(ConfigKey(domainName,"NumGrids"))){
        if(gridNames.empty())
          messageStream << "Number of grids not specified" << std::endl;
        else numGridsConfigured = gridNames.size();
      } else {
        numGridsConfigured = domainConfig.GetValue<int>(ConfigKey(domainName,"NumGrids"));
      }

      if(numGridsConfigured < 0){
        messageStream << "Number of grids cannot be negative." << std::endl;
        return(1);
      } else if (numGridsConfigured == 0){
        messageStream << "Gridless domains are not supported (yet)." << std::endl;
        return(1);
      }

      if(numGridsConfigured > 1){
        messageStream << "Multiple grids detected."
                      << std::endl;
      }

      SetNumGrids(numGridsConfigured);

      bool nativeGeometry = true;
      std::vector<std::string>::iterator gridNameIt = gridNames.begin();
      while(gridNameIt != gridNames.end()){
        int gridIndex = gridNameIt - gridNames.begin();
        std::string &gridName(*gridNameIt++);
        std::string geometryName;
        std::string::size_type geomPos = gridName.find(":");
        if(geomPos != std::string::npos){
          geometryName = gridName.substr(0,geomPos);
        }
        geometryNames.push_back(geometryName);
        //         std::string geomKey(ConfigKey(domainName,gridName+":GeometryName"));
        //         if(domainConfig.IsSet(geomKey)){
        //           nativeGeometry = false;
        //           std::string geometryName(domainConfig.GetValue(geomKey));
        //           std::string geomGridKey(ConfigKey(domainName,gridName+":GridName"));
        //           if(domainConfig.IsSet(geomGridKey)){
        //             geometryName = geometryName + ":" + domainConfig.GetValue(geomGridKey);
        //           }
        //        ownGrid[gridIndex] = true;
        //        domainGrids[gridIndex] = new GridT;
      }
      //       std::sort(geometryNames.begin(),geometryNames.end());
      //       geometryNames.erase(std::unique(geometryNames.begin(),geometryNames.end()),
      //                           geometryNames.end());
      
      // ------ POSTPONE GRID CONFIGURATION ---------------
      //       typename gridvector::iterator gridIt = domainGrids.begin();
      //       while(gridIt != domainGrids.end()){
      //         int iGrid = gridIt - domainGrids.begin();
      //         const std::string &gridName(gridNames[iGrid]);
      //         GridType *domainGridPtr(*gridIt++);
      //         if(domainGridPtr == NULL){
      //           messageStream << "Grid " << iGrid << " failed to configure, invalid storage."
      //                         << std::endl;
      //           return(1);
      //         }
      //         GridType &domainGrid(*domainGridPtr);
      //         if(nativeGeometry){
      //           std::string gridKey(ConfigKey(domainName,gridName));
      //           pcpp::ConfigType gridConfig(GetSubConfig(gridKey,domainConfig));
      //           if(domainGrid->ConfigureGrid(gridConfig,messageStream)){
      //             messageStream << "Grid " << (gridIt - domainGrids.begin())-1
      //                           << " failed to configure." << std::endl;
      //             return(1);
      //           }
      //         }
      //         //         else {
      //         //           if(ConfigureDomainGeometry(domainConfig,geometryNames[iGrid],domainGrid,messageStream)){
      //         //             messageStream << "Failed to configure domain grid " << iGrid << " with geometry "
      //         //                           << geometryNames[iGrid] << "." << std::endl;
      //         //             return(1);
      //         //           }
      //         //         }
      //       }
      
      if(ConfigureDataDictionary(domainConfig,"State",stateDataDictionary,messageStream)){
        messageStream << "Domain state meta failed to configure."
                      << std::endl;
        return(1);
      }
//       std::cout << "Configured state dictionary:" << std::endl
//                 << stateDataDictionary << std::endl;
      if(ConfigureDataDictionary(domainConfig,"Param",paramDataDictionary,messageStream)){
        messageStream << "Domain param meta failed to configure."
                      << std::endl;
        return(1);
      }
//       std::cout << "Configured state dictionary:" << std::endl
//                 << paramDataDictionary << std::endl;
      
      
      // ------ POSTPONE CREATION OF ACTUAL STATE ---------
      //       typename statevector::iterator stateIt = gridStates.begin();
      //       typename statevector::iterator paramIt = gridParams.begin();
      //       while(stateIt != gridStates.end()){
      //         int iGrid = stateIt - gridStates.begin();
      //         StateType *&gridState(*stateIt++);
      //         StateType *&gridParam(*paramIt++);
      //         if(gridState == NULL){
      //           ownState[iGrid] = true;
      //           gridState = new StateT;
      //         }
      //         PopulateGridMeta(gridState->Meta());
      //         if(gridParam == NULL){
      //           gridParam = new StateT;
      //         }
      //         PopulateGridParamMeta(gridParam->Meta());
      //       }
      
      PopulateDomainStateMeta();
      PopulateDomainParamMeta();

      if(domainConfig.IsSet(ConfigKey(domainName,"IO:Fields")))
        ioFieldNames = domainConfig.GetValueVector<std::string>(ConfigKey(domainName,"IO:Fields"));
      
      // Get Domain Boundaries
      bcNames = domainConfig.GetValueVector<std::string>(ConfigKey(domainName,"BCNames"));
      int numBoundaries = bcNames.size();
      //     std::cout << "Number of boundary condition types: " << numBoundaries << std::endl;
      for(int iBoundary = 0;iBoundary < numBoundaries;iBoundary++){

        std::string &bcName(bcNames[iBoundary]);
        //        std::cout << "BCname = " << bcName << std::endl;
        std::string boundaryConfigName(ConfigKey(domainName,bcName));
        //         std::cout << "boundaryType (" << bcName << ") has bc config name: ("
        //                   << boundaryConfigName << ")" << std::endl;
        fixtures::ConfigurationType 
          boundaryConfig(pcpp::util::ExtractConfigParams(domainConfig,boundaryConfigName));

        simulation::domain::boundary::bc::base myBC;
        
        std::string bcTypeName(boundaryConfig.GetValue(ConfigKey(boundaryConfigName,"BCType")));
        //         std::cout << "bcTypeName = " << bcTypeName << std::endl;
        myBC.SetBCName(bcTypeName);
        myBC.SetBCType(boundary::bc::navierstokes::ResolveBCName(myBC.BCName()));
        int bcType = myBC.BCType();
        
        if(bcType != boundary::bc::navierstokes::HOLE){ 
          
          if(ConfigureDataDictionary(boundaryConfig,"State",myBC.StateDictionary(),messageStream)){
            messageStream << "BC state dictionary failed to configure for BC " << bcName 
                          << std::endl;
            return(1);
          }

          if(ConfigureDataDictionary(boundaryConfig,"Param",myBC.ParamDictionary(),messageStream)){
            messageStream << "BC param dictionary failed to configure for BC " << bcName
                          << std::endl;
            return(1);
          }

          if(myBC.ParamDictionary().empty() && myBC.StateDictionary().empty()){
            messageStream << "BC params and fields not configured for " << bcName
                          << std::endl;
            return(1);
          }
          
          // Populate any BC parameter metadata from config, and allocate
          StateType &bcParams(myBC.Params());
          bcParams.SetMetaData(myBC.ParamDictionary());
          bcParams.Create();
          
          // Populate parameter values from config file - any grid-connected state data done later
          std::string bcParamKey(ConfigKey(boundaryConfigName,"Param"));
          if(ConfigStateValues(boundaryConfig,bcParamKey,bcParams,messageStream)){
            messageStream << "BC param values failed to populate from configuration." << std::endl;
            return(1);
          }

        } else {
          messageStream << "Deferring setup of hole region." << std::endl;
        }
        domainBCs.push_back(myBC);
      }
      return(0);
    }
    
    
    int ConfigureGridDomainBoundaries(fixtures::ConfigurationType &inConfig,int iGrid,
                                      std::ostream &inMessageStream)
    {
      
      int numBoundaries = bcNames.size();
      std::ostream &messageStream(inMessageStream);
      
      const std::string &currentGridName(gridNames[iGrid]);
//       std::cout <<  "Processing Grid(" << iGrid << "," 
//                 << currentGridName << ") " << std::endl;
      for(int iBoundary = 0;iBoundary < numBoundaries;iBoundary++){
        
        std::string &bcName(bcNames[iBoundary]);
        //        std::cout << "BCName(" << iBoundary << ") = " << bcName << std::endl;
        BCType &boundaryBC(domainBCs[iBoundary]);
        std::string boundaryConfigKey(ConfigKey(domainName,bcName));
        //        std::cout << "Boundary config key:  " << boundaryConfigKey << std::endl;
        //        std::vector<int> boundaryGridIDs;
        
        // This is a vector of grid regions to which this boundary condition applies
        std::vector<std::string> boundaryGridRegionNames = 
          domainConfig.GetValueVector<std::string>(ConfigKey(boundaryConfigKey,"RegionNames"));
        int numBoundaryGridRegions = boundaryGridRegionNames.size();
        std::vector<std::string>::iterator rnIt = boundaryGridRegionNames.begin();
        //         while(rnIt != boundaryGridRegionNames.end()){
        //           std::cout << "Applies to grid region: " << *rnIt++ << std::endl;
        //         }
        
        boundaryNames = 
          domainConfig.GetValueVector<std::string>(ConfigKey(boundaryConfigKey,"BoundaryNames"));
        rnIt = boundaryNames.begin();
        //         while(rnIt != boundaryNames.end()){
        //           std::cout << "Applies to domain boundary name: " << *rnIt++ << std::endl;
        //         }
        if(boundaryNames.empty() || boundaryNames.size() != boundaryGridRegionNames.size()){
          messageStream << "WARNING: Setting default boundary names." << std::endl;
          boundaryNames.resize(numBoundaryGridRegions);
          for(int iRegion = 0;iRegion < numBoundaryGridRegions;iRegion++)
            boundaryNames[iRegion] = bcName + std::string(":") + boundaryGridRegionNames[iRegion]; 
        }
        // this boundary condition applies to (numBoundaryRegions) grid boundaries
        for(int iRegion = 0;iRegion < numBoundaryGridRegions;iRegion++){
          std::string &boundaryName(boundaryNames[iRegion]);
          std::string gridRegionName(boundaryGridRegionNames[iRegion]);
          //           std::cout << "IRegion(" << iRegion << ") = (" << boundaryName << "," << gridRegionName << ")" 
          //                     << std::endl;
          std::string::size_type x = gridRegionName.find_last_of(":");
          // std::string geometryName;
          std::string regionGridName;
          if(x != std::string::npos){
            regionGridName = gridRegionName.substr(0,x);
            gridRegionName = gridRegionName.substr(x+1);
          }
          if(x == std::string::npos || regionGridName.empty())
            regionGridName = currentGridName;
          std::string regionName = gridRegionName;
          int gridIndex = 0;
          //           std::cout << "RegionName: " << regionName << std::endl;
          //           std::cout << "regionGridName: " << regionGridName << std::endl;
          if(regionGridName != currentGridName){
            if(!regionGridName.empty()){ // make sure *some* grid region will match
              gridIndex = GetGridIndex(regionGridName);
              if(gridIndex < 0){
                messageStream << "Failed to find grid id for \""
                              << regionGridName << "\" in BC resolve."
                              << std::endl;
                return(1);
              }
            }
            continue;
          }
            
          gridIndex = iGrid;
            
          std::vector<BoundaryType> &domainGridBoundaries(GridBoundaries(gridIndex));
          std::string &bcName(boundaryBC.BCName());
          //          std::cout << "BCTypeName: " << bcName << std::endl;
          //  domainConfig.GetValue<std::string>(ConfigKey(boundaryConfigKey,"Type")));
          //  myBoundary.SetBCName(domainConfig.GetValue<std::string>(ConfigKey(boundaryConfigKey,"Type")));
          GridType &boundaryGrid(Grid(gridIndex));
          //           std::cout << "gridIndex = " << gridIndex << std::endl;
          //           std::cout << "Finding grid region for regionName = " << regionName
          //                     << std::endl; 
          int gridRegionIndex = boundaryGrid.GetRegionIndex(regionName);
          //          std::cout << "gridRegionIndex = " << gridRegionIndex << std::endl;
          if(gridRegionIndex < 0){
            messageStream << "Failed to find grid region id for \""
                          << gridNames[gridIndex] << ":" << regionName 
                          << "\" in BC resolve." << std::endl;
            //            std::cout << "Did not find region index, aborting." << std::endl;
            return(1);
          }
            
          BoundaryType thisBoundary;
          thisBoundary.SetGridRegionID(gridRegionIndex);
          thisBoundary.SetGrid(boundaryGrid);
          thisBoundary.SetBC(boundaryBC);
          thisBoundary.SetName(boundaryName);
          messageStream << "Found boundary " << boundaryName << " with BC " 
                        << bcName << " for (grid:region) = ("
                        << gridNames[gridIndex] << ":" << regionName
                        << ")" << std::endl;
          //          std::cout << "Adding boundary!" << std::endl;
          domainGridBoundaries.push_back(thisBoundary);
        }
      }      
      return(0);
    };

    int CreateData(std::ostream &messageStream)
    {
      int numGrids = domainGrids.size();
      int numGridStates = gridStates.size();
      int numGridParams = gridParams.size();

      if(numGrids != numGridStates)
        return(1);
      if(numGridStates != numGridParams)
        return(1);

      std::string stateFields;
      if(domainConfig.IsSet("Advancer:Fields")){
        stateFields = domainConfig.GetValue(ConfigKey(domainName,"Advancer:Fields"));
        messageStream << "Setting state fields to (" << stateFields << ")" << std::endl;
      } else {
        messageStream << "WARNING: No fields set up as State Fields." << std::endl;
      }
      typename gridvector::iterator gridIt   = domainGrids.begin();
      typename statevector::iterator stateIt = gridStates.begin();
      typename statevector::iterator targIt  = gridTargets.begin();
      typename statevector::iterator paramIt = gridParams.begin();
      std::vector<std::string>::iterator gridNameIt = gridNames.begin();

      messageStream << "Creating state and parameter data for " << domainGrids.size() << " grids."
                    << std::endl;

      while(gridIt != domainGrids.end()){

        GridType  &thisGrid(**gridIt++);
        StateType &gridState(**stateIt++);
        //        StateType &gridTarget(**targIt++);
        StateType &gridParams(**paramIt++);
        std::string &gridName(*gridNameIt++);

        const std::vector<size_t> &bufferSizes(thisGrid.BufferSizes());
        const pcpp::IndexIntervalType &partitionInterval(thisGrid.PartitionInterval());
        const pcpp::IndexIntervalType &partitionBufferInterval(thisGrid.PartitionBufferInterval());  
        size_t numPointsBuffer = thisGrid.BufferSize();


        messageStream << "Creating buffers for " << gridName
                      << " with bufferSizes(";
        pcpp::io::DumpContents(messageStream,bufferSizes,",");
        messageStream << ")" << std::endl;

        gridState.SetMetaData(stateDataDictionary);
        gridParams.SetMetaData(paramDataDictionary);

        gridState.Create(numPointsBuffer,0);
        if(!stateFields.empty()) 
          gridState.SetStateFields(stateFields);
        if(*targIt){
          (*targIt)->CopyStateFields(gridState);
        }
        gridParams.Create(numPointsBuffer,0);

        std::string paramKeyRoot(ConfigKey(domainName,"Param"));
        messageStream << "Configuring parameters with KeyRoot(" 
                      << paramKeyRoot << ")." << std::endl;
        if(ConfigStateValues(domainConfig,paramKeyRoot,gridParams,messageStream)){
          messageStream << "ERROR! Domain parameters failed to populate." << std::endl;
          return(1);
        }
      }
      return(0);
    };

    int CreateData(int gridIndex,std::ostream &messageStream)
    {
      if(gridIndex > (domainGrids.size()-1))
        return(1);

      std::string stateFields(domainConfig.GetValue(ConfigKey(domainName,"Advancer:Fields")));
      messageStream << "Setting state fields to (" << stateFields << ")" << std::endl;


      messageStream << "Creating state and parameter data for grid(" << gridIndex << ","
                    << gridNames[gridIndex] << ")" << std::endl;

      
      GridType  &thisGrid(Grid(gridIndex));
      StateType &gridState(State(gridIndex));
      //      StateType &gridTarget(Target(gridIndex));
      StateType &gridParams(Param(gridIndex));
      std::string &gridName(gridNames[gridIndex]);
      
      const std::vector<size_t> &bufferSizes(thisGrid.BufferSizes());
      const pcpp::IndexIntervalType &partitionInterval(thisGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(thisGrid.PartitionBufferInterval());  
      size_t numPointsBuffer = thisGrid.BufferSize();
      
      
      messageStream << "Creating buffers for " << gridName
                    << " with bufferSizes(";
      pcpp::io::DumpContents(messageStream,bufferSizes,",");
      messageStream << ")" << std::endl;

      gridState.SetMetaData(stateDataDictionary);
      gridParams.SetMetaData(paramDataDictionary);
      
      gridState.Create(numPointsBuffer,0);
      gridState.SetStateFields(stateFields);
      if(gridTargets[gridIndex] == NULL){
        gridTargets[gridIndex] = new StateType;
      }
      gridTargets[gridIndex]->CopyStateData(gridState);

      gridParams.Create(numPointsBuffer,0);
      
      std::string paramKeyRoot(ConfigKey(domainName,"Param"));
      if(ConfigStateValues(domainConfig,paramKeyRoot,gridParams,messageStream)){
        messageStream << "ERROR! Domain parameters failed to populate." << std::endl;
        return(1);
      }

      return(0);

    };

    int GetGridIndex(const std::string &gridName) const
    {
      int gridIndex = 0;
      std::vector<std::string>::const_iterator nameIt = gridNames.begin();
      while(nameIt != gridNames.end()){
        if(*nameIt++ == gridName){
          return(gridIndex);
        }
        gridIndex++;
      }
      return(-1);
    };


    int NumberOfGrids(){ return(numGrids); };
    int NumberOfLocalGrids(){ return(partitionInfo.gridIndexOfLocalGrids.size()); };
    std::vector<int> &LocalGridIndices()
    { 
      return(partitionInfo.gridIndexOfLocalGrids); 
    };
    int LocalGridIndex(int localGridIndex)
    { 
      return(partitionInfo.gridIndexOfLocalGrids[localGridIndex]); 
    };

    std::string Report()
    {
      std::ostringstream Ostr;
      Ostr << "+++Domain data information+++ " << std::endl
           << "State data dictionary:        " << std::endl
           << stateDataDictionary.Report()     << std::endl
           << "Parameter data dictionary:    " << std::endl
           << paramDataDictionary.Report()     << std::endl;
      if(!domainBCs.empty()){
        Ostr << "Boundary Conditions:          " << std::endl;
        std::vector<BCType>::iterator bcIt = domainBCs.begin();
        int bcIndex = 0;
        while(bcIt != domainBCs.end()){
          BCType &myBC(*bcIt++);
          Ostr << "BC(" << bcNames[bcIndex++] << ") = " 
               << myBC.BCName() << std::endl;
          if(!myBC.StateDictionary().empty()){
            Ostr << "  Boundary state data: " << std::endl
                 << myBC.StateDictionary().Report() << std::endl;
          } else {
            Ostr << "  No boundary state data. " << std::endl;
          }
          if(!myBC.ParamDictionary().empty()){
            Ostr << "  Boundary params: " << std::endl
                 << myBC.ParamDictionary().Report() << std::endl;
          } else {
            Ostr << "  No boundary parameters. " << std::endl;
          }
        }
      } else {
        Ostr << "No boundary conditions configured."      << std::endl;
      }
      Ostr << "Number of grids: " << numGrids  << std::endl
           << "... Grids ..." << std::endl;
      for(int iGrid = 0;iGrid < numGrids;iGrid++){
        Ostr << "Grid(" << iGrid+1 << "):" << std::endl
             << "Grid Name: " << gridNames[iGrid] << std::endl
             << "Grid Geometry: " << geometryNames[iGrid] << std::endl;
        if(domainGrids[iGrid] != NULL){
          Ostr << "Grid Configuration: " << std::endl
               << domainGrids[iGrid]->ReportConfiguration() << std::endl;
          std::vector<size_t> &gridSizes(domainGrids[iGrid]->GridSizes());
          std::vector<size_t> &bufferSizes(domainGrids[iGrid]->BufferSizes());
          pcpp::IndexIntervalType &partitionInterval(domainGrids[iGrid]->PartitionInterval());
          pcpp::IndexIntervalType &partitionBufferInterval(domainGrids[iGrid]->PartitionBufferInterval());
          size_t numPointsBuffer = domainGrids[iGrid]->BufferSize();
          Ostr << "Grid Sizes: (";
          pcpp::io::DumpContents(Ostr,gridSizes,",");
          Ostr << ")" << std::endl
               << "Partition Interval: ";
          partitionInterval.PrettyPrint(Ostr);
          Ostr << std::endl
               << "Buffer Sizes: (";
          pcpp::io::DumpContents(Ostr,bufferSizes,",");
          Ostr << ")" << std::endl
               << "Partition Buffer Interval: ";
          partitionBufferInterval.PrettyPrint(Ostr);
          Ostr << std::endl
               << "Buffer size: " << numPointsBuffer << std::endl;
        } 
      }
      return(Ostr.str());
    };
    DataDictionaryType &StateDictionary() { return(stateDataDictionary); };
    DataDictionaryType &ParamDictionary() { return(paramDataDictionary); };

    void PopulateGridMeta(DataDictionaryType &gridMeta)
    {     
      typename DataDictionaryType::iterator metaIt = stateDataDictionary.begin();
      while(metaIt != stateDataDictionary.end()){
        typename StateType::MetaDataType &metaData(*metaIt++);
        if(metaData.loc == 'g' || metaData.loc == 'm' ||
           metaData.loc == 'c' || metaData.loc == 'n')
          gridMeta.AddField(metaData.name,metaData);
      }
    };

    void PopulateGridParamMeta(DataDictionaryType &gridMeta)
    {     
      typename DataDictionaryType::iterator metaIt = paramDataDictionary.begin();
      while(metaIt != paramDataDictionary.end()){
        typename StateType::MetaDataType &metaData(*metaIt++);
        if(metaData.loc == 'g' || metaData.loc == 'm' ||
           metaData.loc == 'c' || metaData.loc == 'n')
          gridMeta.AddField(metaData.name,metaData);
      }
    };

    void PopulateDomainStateMeta()
    {     
      typename DataDictionaryType::iterator metaIt = stateDataDictionary.begin();
      while(metaIt != stateDataDictionary.end()){
        typename StateType::MetaDataType &metaData(*metaIt++);
        if(metaData.loc == 'g' || metaData.loc == 'd' ||
           metaData.loc == 's')
          domainState.Meta().AddField(metaData.name,metaData);
        //          gridMeta.AddField(metaData.name,metaData);
      }
    };
    
    void PopulateDomainParamMeta()
    {     
      typename DataDictionaryType::iterator metaIt = paramDataDictionary.begin();
      while(metaIt != paramDataDictionary.end()){
        typename StateType::MetaDataType &metaData(*metaIt++);
        if(metaData.loc == 'g' || metaData.loc == 'd' ||
           metaData.loc == 's')
          domainParams.Meta().AddField(metaData.name,metaData);
        //          gridMeta.AddField(metaData.name,metaData);
      }
    };

    fixtures::ConfigurationType &Config(){return(domainConfig);};
    std::vector<std::string> &GeometryNames() {return(geometryNames);};
    std::vector<std::string> &GridNames() {return(gridNames); };
    std::string GridName(int iGrid) {return(gridNames[iGrid]); };
    std::vector<std::string> IOFieldNames() { return(ioFieldNames); };
    void SetSpatialOrder(int inOrder){ spatialOrder = inOrder;};
    int SpatialOrder(){return(spatialOrder);};
    void SetTime(double inTime){ domainTime = inTime; };
    double Time(){return(domainTime); };
    OperatorType &Operator(){return(domainOperator); };
    const OperatorType &Operator() const {return(domainOperator); };
    void SetMessageStream(std::ostream &inStream){messageStreamPtr = &inStream;};
    std::vector<BCType> &BCs() { return(domainBCs); };
    const std::vector<BCType> &BCs() const { return(domainBCs); };
    RHSType &RHS() { return(*rhsPtr); };
    int SetRHS(RHSType &inRHS) 
    { 
      if(numGrids <= 0){
        return(1);
      }
      if(SetRHS(0,inRHS))
        return(1);
      rhsPtr = gridRHSs[0]; 
      return(0);
    };
    int SetRHS(int gridIndex,RHSType &inRHS) 
    { 
      if(partitionInfo.gridIndexOfLocalGrids.size() != numGrids){
        if(PartitionDomain(0)){
          return(1);
        }
      }
      if(gridRHSs.size() != numGrids)
        gridRHSs.resize(numGrids,NULL);
      if(ownRHS.size() != numGrids)
        ownRHS.resize(numGrids,false);
      if(ownRHS[gridIndex] && (gridRHSs[gridIndex] != NULL)){
        delete gridRHSs[gridIndex];
        ownRHS[gridIndex] = false;
      }
      gridRHSs[gridIndex] = &inRHS;
      return(0);
    };
    fixtures::CommunicatorType &Communicator() { return(domainCommunicator); };
    fixtures::CommunicatorType &GridCommunicator(int gridIndex) 
    {
      return(gridCommunicators[gridIndex]); 
      //      return(domainGrids[gridIndex]->Communicator()); 
    };
    std::vector<fixtures::CommunicatorType> &GridCommunicators() {return(gridCommunicators); };
    fixtures::CommunicatorType &SetCommunicator(fixtures::CommunicatorType &inCommunicator) 
    { return(domainCommunicator = inCommunicator); };
    int PartitionDomain(int inOptions = 0)
    { 
      int returnCode = 0;
      if(inOptions == 0){ // Trivial "partitioning" where all grids are on all procs
        int numLocalGrids = domainGrids.size();
        partitionInfo.gridIndexOfLocalGrids.resize(numLocalGrids);
        for(int iGrid = 0;iGrid < numLocalGrids;iGrid++){
          partitionInfo.gridIndexOfLocalGrids[iGrid] = iGrid;
        }
        //        gridCommunicators.resize(numLocalGrids,domainCommunicator);
      } else if(inOptions == 1){ // Auto-split grids according to size
        int numDomainProcs = 1;
        if(domainCommunicator.Good()){
          numDomainProcs = domainCommunicator.NProc();
          domainRank = domainCommunicator.Rank();
        }
        if(numDomainProcs == 1 || numGrids > numDomainProcs)
          return(PartitionDomain(0));
        size_t numPointsTotal = 0;
        std::vector<size_t> numPointsGrid(numGrids,0);
        std::vector<int>  numProcsGrid(numGrids,0);
        std::vector<size_t> numPointsPerProc(numGrids,0);
        for(int iGrid = 0;iGrid < numGrids;iGrid++){
          GridType &currentGrid(*domainGrids[iGrid]);
          const std::vector<size_t> &gridSizes(currentGrid.GridSizes());
          int numDim = gridSizes.size();
          for(int iDim = 0;iDim < numDim;iDim++){
            numPointsGrid[iGrid] += gridSizes[iDim];
          }
          numPointsTotal += numPointsGrid[iGrid];
        }
        size_t targetPointsPerProc = numPointsTotal/numDomainProcs;
        int totalGridProcs = 0;
        for(int iGrid = 0;iGrid < numGrids;iGrid++){
          numProcsGrid[iGrid] = int(numPointsGrid[iGrid]*numDomainProcs/float(numPointsTotal));
          numProcsGrid[iGrid] = std::max(1,numProcsGrid[iGrid]);
          totalGridProcs += numProcsGrid[iGrid];
          numPointsPerProc[iGrid] = numPointsGrid[iGrid]/numProcsGrid[iGrid];
        }
        if(totalGridProcs > numDomainProcs){
          int numToScrape = totalGridProcs - numDomainProcs;
          std::vector<unsigned int> numProcsOrder;
          ix::util::SortPermutation(numProcsGrid,numProcsOrder);
          std::reverse(numProcsOrder.begin(),numProcsOrder.end());
          int orderIndex = 0;
          for(int iScrape = numToScrape;iScrape > 0;iScrape--){
            int gridIndexToScrape = numProcsOrder[orderIndex];
            numProcsGrid[gridIndexToScrape]--;
            numPointsPerProc[gridIndexToScrape] = 
              numPointsGrid[gridIndexToScrape]/numProcsGrid[gridIndexToScrape];
          }
        } else {
          int numToDoleOut = numDomainProcs - totalGridProcs;
        } 
        //        if(totalGridProcs == numDomainProcs)
          //          return(PartitionDomain(numProcsGrid));

      } else if(inOptions == 3){  // Fully user prescribed domain splitting

        if(!geometryInfoPtr)
          return(1);
       //        if((splitMap.size() != numGrids) || slitMap.empty())
        //           return(1);
        
        int numDomainProcs = 1;
        if(domainCommunicator.Good()){
          numDomainProcs = domainCommunicator.NProc();
          domainRank = domainCommunicator.Rank();
        }

        int numProcsTotal = 0;
        size_t numPointsTotal = 0;
        std::vector<size_t> numPointsGrid(numGrids,1);
        std::vector<int>  numProcsGrid(numGrids,1);
        std::vector<size_t> numPointsPerProc(numGrids,0);
        for(int iGrid = 0;iGrid < numGrids;iGrid++){
          GridType &currentGrid(*domainGrids[iGrid]);
          const std::vector<size_t> &gridSizes(currentGrid.GridSizes());
          int numDim = gridSizes.size();
          const std::vector<int> &gridSplitSizes(currentGrid.DecompSizes());
          if(gridSplitSizes.size() != numDim)
            return(1);          
          for(int iDim = 0;iDim < numDim;iDim++){
            numProcsGrid[iGrid]  *= gridSplitSizes[iDim];
            numPointsGrid[iGrid] *= gridSizes[iDim];
          }
          if(numPointsGrid[iGrid] == 0)
            return(1);
          if(!(numProcsGrid[iGrid] > 0)) 
            return(1);
          numPointsTotal += numPointsGrid[iGrid];
          numProcsTotal  += numProcsGrid[iGrid];
          numPointsPerProc[iGrid] = numPointsGrid[iGrid]/numProcsGrid[iGrid];
          if(numPointsPerProc[iGrid] == 0 || numPointsPerProc[iGrid] > numPointsGrid[iGrid])
            return(1);
        }
        if(numProcsTotal != numDomainProcs)
          return(1);
        
        // Now split communicator into grid-specific ones
        //        gridCommunicators.resize(numGrids,domainCommunicator);
        int startRank = 0;
        for(int iGrid = 0;iGrid < numGrids;iGrid++){
          int endRank = startRank + numProcsGrid[iGrid];
          int color = ((domainRank >= startRank) && (domainRank < endRank));
          domainCommunicator.Split(color,domainRank,gridCommunicators[iGrid]);
          startRank = endRank;
        }
      } else {
        returnCode = 1;
      }
      
      return(returnCode);
    };

    int InitializeRHS(int inOptions)
    {
      if(gridRHSs.size() != numGrids){
        gridRHSs.resize(numGrids,NULL);
      }
      if(ownRHS.size() != numGrids){
        ownRHS.resize(numGrids,false);
      }
      std::vector<int>::iterator gridIdIt = partitionInfo.gridIndexOfLocalGrids.begin();
      while(gridIdIt != partitionInfo.gridIndexOfLocalGrids.end()){
        int gridIndex = *gridIdIt++;
        gridRHSs[gridIndex] = new RHSType;
        ownRHS[gridIndex] = true;
        RHSType &gridRHS(*gridRHSs[gridIndex]);
        
      }
    };
    
    int Step() { return(stepNumber); };
    void SetStep(int inStep) { stepNumber = inStep; };

    // Must be called on all threads
    double TimeStep(double *outCFL = NULL)
    {
      double timeStep = std::numeric_limits<double>::max();
      double cfl = 1.0;
      //      if(timeSteppingMode == CONSTANT_DT){
      if(true){
        std::vector<int>::iterator gridIdIt = partitionInfo.gridIndexOfLocalGrids.begin();
        while(gridIdIt != partitionInfo.gridIndexOfLocalGrids.end()){
          int gridIndex = *gridIdIt++;
          double gridDT = gridRHSs[gridIndex]->TimeStep();
          if(gridDT < timeStep) timeStep = gridDT;
        }
      }
      return(timeStep);
    };

    // Must be called on all threads
    void ComputeDV(int threadId = 0)
    {
      std::vector<int> &gridIndexOfLocalGrids(partitionInfo.gridIndexOfLocalGrids);
      std::vector<int>::iterator gridIdIt = gridIndexOfLocalGrids.begin();
      while(gridIdIt != gridIndexOfLocalGrids.end()){
        int gridIndex = *gridIdIt++;
        gridRHSs[gridIndex]->ComputeDV(threadId);
      }
    };

    int SetNumThreads(int numThreads)
    {
      std::vector<int> &gridIndexOfLocalGrids(partitionInfo.gridIndexOfLocalGrids);
      std::vector<int>::iterator gridIdIt = gridIndexOfLocalGrids.begin();
      while(gridIdIt != gridIndexOfLocalGrids.end()){
        int gridIndex = *gridIdIt++;
        domainGrids[gridIndex]->SetNumThreads(numThreads);
        gridRHSs[gridIndex]->InitThreadIntervals();
      }
    };

    int SetDomainBCs()
    {
      std::vector<int> &gridIndexOfLocalGrids(partitionInfo.gridIndexOfLocalGrids);
      std::vector<int>::iterator gridIdIt = gridIndexOfLocalGrids.begin();
      while(gridIdIt != gridIndexOfLocalGrids.end()){
        int gridIndex = *gridIdIt++;
        gridRHSs[gridIndex]->SetDomainBCs(DomainBoundary(gridIndex),BCs());
      }
      return(0);
    };

    int InitThreads()
    {
      std::vector<int> &gridIndexOfLocalGrids(partitionInfo.gridIndexOfLocalGrids);
      std::vector<int>::iterator gridIdIt = gridIndexOfLocalGrids.begin();
      while(gridIdIt != gridIndexOfLocalGrids.end()){
        int gridIndex = *gridIdIt++;
        gridRHSs[gridIndex]->InitThreadIntervals();
      }
      return(0);
    };
    
    int FinalizeAdvance(double newTime){
      int advancementMode = 0;
      if(advancementMode == 0){ // just barrier until all grids are done
        // DOH!  this caused failure because MPI not initialized - added a
        // check to see if the domain communicator is "good" before 
        // invoking barrier.
        if(domainCommunicator.Good())
          domainCommunicator.Barrier();
        // could/should validate here that all grids have the same
        // time - at least for this mode.
        SetTime(newTime);
        return(0);
      } else {
        return(1);
      }
      return(0);
    };
 
    void SetGeometryInfo(geometry::info &geometryInfo)
    { geometryInfoPtr = &geometryInfo; };
    
    void SetScalarNames(const std::string &inNames){ scalarNames = inNames; };
    std::vector<std::string> &ScalarNames() { return(scalarNames); };
    const std::vector<std::string> &ScalarNames() const { return(scalarNames); };

  protected:

    geometry::info                          *geometryInfoPtr;
    pcpp::ParallelGlobalType                *globalPtr;
    std::ostream                            *messageStreamPtr;
    std::string                             domainName;
    OperatorType                            domainOperator;    
    StateType                               domainState;
    StateType                               domainParams;
    gridvector                              domainGrids;
    statevector                             gridStates;
    statevector                             gridTargets;
    statevector                             gridParams;
    rhsvector                               gridRHSs;
    int                                     spatialOrder;
    int                                     stepNumber;
    double                                  domainTime;
    std::vector<bool>                       ownRHS;
    std::vector<bool>                       ownGrid;
    std::vector<bool>                       ownState;
    fixtures::ConfigurationType             domainConfig;
    DataDictionaryType                      stateDataDictionary;
    DataDictionaryType                      paramDataDictionary;
    std::vector<std::string>                gridNames;
    std::vector<std::string>                geometryNames;
    std::vector<std::string>                ioFieldNames;
    std::vector<std::vector<BoundaryType> > domainBoundaries;
    std::vector<std::string>                bcNames;
    std::vector<std::string>                boundaryNames;
    std::vector<std::string>                scalarNames;
    std::vector<BCType>                     domainBCs;
    RHSType                                 *rhsPtr;
    RHSType                                 dummyRHS;
    fixtures::CommunicatorType domainCommunicator;
    std::vector<fixtures::CommunicatorType> gridCommunicators;
    int domainRank;
    int                                 numGrids;
    partitioninfo                       partitionInfo;

  };
  
  
  template<typename DomainType>
  class Initializer {
  protected:
    fixtures::ConfigurationType initConfig;
  public:
    virtual int Configure(fixtures::ConfigurationType &inConfig) {return(0);};
    virtual int InitializeDomain(DomainType &inDomain) {return(0);};
  };
  
  // namespace initializer {
  //   class base {
  //   protected:
  //     ConfigurationType initConfig;
  //   public:
  //     virtual int Configure(ConfigurationType &inConfig)  { return (0); };
  //     template<typename GridT,typename StateT>
  //     virtual int InitializeDomain(domain::base<GridT,StateT> &inDomain){ return (0); };
  //   };
    
  // };
};
};
#endif