PC2Configuration.C

Go to the documentation of this file.

#include "PlasCom2.H"
#include "PCPPReport.H"
#include "PCPPUtil.H"
#include "PCPPIntervalUtils.H"
#include "PC2DefaultConfiguration.H"
#include "PC2Util.H"
#include "Report.H"
#include "GridGenerator.H"

namespace plascom2 {

  namespace util {

    int ConfigureVirtualNode(const fixtures::ConfigurationType &inConfig,
                             fixtures::CommunicatorType &inCommunicator,
                             pcpp::VirtualNodeType &virtualNode,
                             std::ostream &messageStream)
    {
      int iopn = 0;
      int ppn  = 0;
      if(inConfig.IsSet("VirtualNode:NumIO"))
        iopn = inConfig.GetValue<int>("VirtualNode:NumIO");
      if(inConfig.IsSet("VirtualNode:NumProc"))
        ppn  = inConfig.GetValue<int>("VirtualNode:NumProc");
      int myRank = inCommunicator.Rank();
      int numProc = inCommunicator.Size();
      if(ppn <= 0) ppn = 1;
      if(iopn < 0) iopn = 0;
      if(ppn == 1){
        virtualNode.procsPerNode = inCommunicator.Size();
      } else {
        virtualNode.procsPerNode = ppn;
      }
      int numNodes = numProc/virtualNode.procsPerNode;
      int nodeColor = myRank/virtualNode.procsPerNode;
      virtualNode.numNodes = numNodes;
      inCommunicator.Split(nodeColor,myRank,virtualNode.nodeCommunicator);
      virtualNode.nodeRank = virtualNode.nodeCommunicator.Rank();
      virtualNode.performIO      = 1;
      virtualNode.performCompute = 1;
      if(iopn > 0) {
        virtualNode.ioPerNode = iopn;
        if(virtualNode.nodeRank > 0)
          virtualNode.performIO = 0;
        else
          virtualNode.performCompute = 0;
      }
      if(virtualNode.performIO > 0 && 
         virtualNode.performCompute > 0){
        virtualNode.ioCommunicator      = inCommunicator;
        virtualNode.computeCommunicator = inCommunicator;
        virtualNode.ioRank              = myRank;
        virtualNode.computeRank         = myRank;
      } else if(virtualNode.performCompute > 0){
        inCommunicator.Split(virtualNode.performCompute,myRank,virtualNode.computeCommunicator);
        virtualNode.ioRank      = -1;
        virtualNode.computeRank = virtualNode.computeCommunicator.Rank();
      } else {
        inCommunicator.Split(virtualNode.performCompute,myRank,virtualNode.ioCommunicator);
        virtualNode.computeRank = -1;
        virtualNode.ioRank      = virtualNode.ioCommunicator.Rank();
      }
      return(0);
    };

    int SplitDomains(const fixtures::ConfigurationType &inConfig,
                     fixtures::CommunicatorType &inCommunicator, 
                     application_t::DomainInfoType &domainPartitionInfo,
                     std::vector<application_t::DomainBaseType> &localDomains,
                     std::ostream &messageStream)
    {
      int numDomains = domainPartitionInfo.numDomains;
      int numProc     = inCommunicator.Size();
      if(numDomains == 0) { 
        messageStream << "No domains to split." << std::endl;
        return(1); 
      } else if(numDomains == 1 || numDomains > numProc) {
        domainPartitionInfo.numLocalDomains = numDomains;
        domainPartitionInfo.globalDomainIndex.resize(numDomains,0);
        localDomains.resize(numDomains);
        for(int i = 0;i < numDomains;i++){
          localDomains[i].SetCommunicator(inCommunicator);
          domainPartitionInfo.globalDomainIndex[i] = i;
        }
        messageStream << (numProc > 1 ? "Splitting " : "Assigning ")
                      << numDomains << " domain" 
                      << (numDomains > 1 ? "s " : " ")
                      << (numProc > 1 ? "across " : "to ")
                      << numProc << " processor" 
                      << (numProc > 1 ? "s." : ".")
                      << std::endl;
        return(0);
      }
      if(domainPartitionInfo.domainNames.size() != numDomains){
        messageStream << "Missing domain names." << std::endl;
        return(1);
      }
      int myRank = inCommunicator.Rank();
      int procsPerDomain = numProc/numDomains;
      if(procsPerDomain <= 0){
        messageStream << "Processors per domain is 0!?" << std::endl;
        return(1);
      }
      int leftOvers = numProc%numDomains;
      if(leftOvers > 0)
        messageStream << "Domain splitting is lopsided with " << leftOvers 
                      << " leftover processes." << std::endl;
      int color = -1;
      int endingRank = 0;
      for(int i = 0;i < numDomains;i++){
        int startingRank = endingRank;
        endingRank   += procsPerDomain;
        if(i < leftOvers) endingRank++;
        if(myRank <= startingRank && myRank < endingRank){
          color = i;
          domainPartitionInfo.numLocalDomains = 1;
          domainPartitionInfo.globalDomainIndex.resize(1,i);
          localDomains.resize(1);
          messageStream << "Rank " << myRank << " taking domain "
                        << i << ": " << domainPartitionInfo.domainNames[i]
                        << std::endl;
        }
      }
      inCommunicator.Split(color,myRank,localDomains[0].Communicator());
      fixtures::CommunicatorType &domainCommunicator(localDomains[0].Communicator());
      int domainRank = domainCommunicator.Rank();
      int numDomainProcs = domainCommunicator.Size();
      messageStream << "DomainRank (global,domain/domain size):"  << myRank 
                    << "," << domainRank << "/" << numDomainProcs << std::endl;
      return(0);
    }
    
    int ConfigurePlatformInfo(const fixtures::ConfigurationType &inConfig,
                              pcpp::PlatformInfoType &platformInfo)
    {
      return(0);
    }
  }
  
  template<typename GridType>
  int GenerateGrid(GridType &inGrid,const std::string &gridName,
                   const fixtures::ConfigurationType &inConfig,
                   std::ostream &messageStream)
  {
    const std::string gridGenMode(inConfig.GetValue(gridName+":GeneratorMode"));

    if(gridGenMode.empty()){
      messageStream << "GenerateGrid:Error: No generation mode specified." << std::endl;
      return(1);
    }
    
    if(gridGenMode == "Cylinder"){
      int gridDimension = inGrid.Dimension();
      if(gridDimension == 2) {
        messageStream << "GenerateGrid: Generating grid with Cylinder2D." << std::endl;
        std::vector<double> gridGenParameters(inConfig.GetValueVector<double>(gridName+":GeneratorParameters"));
        std::vector<int>    gridGenOptions(inConfig.GetValueVector<int>(gridName+":GeneratorOptions"));
        if(gridGenParameters.empty()){
          gridGenParameters.resize(2,.5);
          gridGenParameters[1] = 1.5;
        }
        inGrid.SetGridGenerationParameters(gridGenParameters);
        inGrid.SetGridGenerationOptions(gridGenOptions);
        return(inGrid.GenerateGrid(gridgen::Cylinder2D,messageStream));
      } else {
        messageStream << "GenerateGrid:Error: 3D Cylinders not supported yet." << std::endl;
        return(1);
      }
    } else {
      messageStream << "GenerateGrid:Error: Unknown grid generation mode ("
                    << gridGenMode << ")" << std::endl;
      return(1);
    }
    return(1);
  };

  int application::ConfigureApplication()

  {
    FunctionEntry("Configure");
    std::ostringstream messageStream;

    int myRank = 0;
    if(globalCommunicator.Good())
      myRank = globalCommunicator.Rank();
    
    // Populate the configuration object from the config file, or restart file
    std::string restartFileName; // temporarily set the restartFileName here
    std::vector<std::string> pc2Arguments(CommandLine().GetArgs());

    if(!pc2Arguments.empty()){
      restartFileName = pc2Arguments[0];
      SetRestart();
      fixtures::ConfigurationType::ParamType param;
      param.first="RestartFileName";
      param.second=restartFileName;
      appConfig.push_back(param);
      if(pc2Arguments.size() > 1){
        param.first="TargetFileName";
        param.second=pc2Arguments[1];
        appConfig.push_back(param);
      }
    }

    std::string configFileName(CommandLine().GetOption('c'));
    
    if(!configFileName.empty()){
      StdOut(std::string("Configuration File Name: ")+configFileName+std::string("\n"),2);
      if(pcpp::util::ProcessConfigurationFile(configFileName,appConfig,globalCommunicator,messageStream)){
        ErrOut(messageStream.str()+std::string("\n"));
        ErrOut(std::string("Processing configuration file ("+configFileName+std::string(") failed.\n")));
        return(2);
      } else {
        //         StdOut(messageStream.str(),2);
        pcpp::io::RenewStream(messageStream);
        //   if(Configure())
        //     return(1);
      }
    } else if (!restartFileName.empty()) { // no configuration file specified
      // get config from restart file (do it now)
      if(pcpp::util::ProcessConfigurationRestart(restartFileName,appConfig,globalCommunicator,messageStream)){
        ErrOut(messageStream.str()+std::string("\n"));
        ErrOut(std::string("Processing configuration from restart file (")+
               restartFileName+std::string(") failed.\n"));
        return(2);
      }
    } else {
      // Set the defaults 
      std::istringstream Istr(PC2DefaultConfiguration);
      Istr >> appConfig;
    } 
    
    // Add any build information 
    ConfigureBuildInformation(appConfig);

    // Should stop if configuration not possible
    std::string platformInfoString;
    platformInfo.hostName.assign(globalCommunicator.Hostname());
    if(appConfig.GetFlagValue("System:Discovery")){
      if(globalCommunicator.Rank() == 0) {
        platformInfo = pcpp::PlatformInfo(pcpp::SystemInfo());
        platformInfoString = platformInfo.Pack();
      }
      globalCommunicator.BroadCast(platformInfoString);
      platformInfo.UnPack(platformInfoString);
    } else {
      // Get platform info from configuration 
      //      if(util::ConfigurePlatformInfo(appConfig,platformInfo))
        
    }

    // Report the platform information (if verbosity high enough)
    std::ostringstream processInfoOut;
    processInfoOut << "Platform Information (" 
                   << globalCommunicator.Rank() << "/"
                   << globalCommunicator.Size() << "):"  << std::endl
                   << "--------------------------------" << std::endl
                   << platformInfo                       << std::endl;
    StdOut(processInfoOut.str(),2);


    // Convert the legacy PlasComCM configuration to PlasCom2 config
    if(appConfig.IsSet("FLUID_MODEL")){
      if(util::ConvertLegacyConfiguration(appConfig,globalCommunicator,messageStream)){
        ErrOut("Failed to convert legacy configuration. Messages:\n");
        ErrOut(messageStream.str());
        return(1);
      }
      if(myRank == 0){
        std::ofstream Ouf;
        Ouf.open("plascom2.config");
        Ouf << appConfig << std::endl;
        Ouf.close();
      }
    }
    
    globalCommunicator.Barrier();
    
    pcpp::report::Configuration(messageStream,appConfig);
    StdOut(messageStream.str(),3);
    pcpp::io::RenewStream(messageStream);

    globalCommunicator.Barrier();

    // If there is a simulation to do, then configure it
    if(appConfig.IsSet("Domain:Names")){
      int simulationConfigCode = ConfigureSimulation();
      if(simulationConfigCode){
        std::ostringstream errStream;
        errStream << "Simulation failed to configure with error code (" 
                  << simulationConfigCode << ")." << std::endl;
        ErrOut(errStream.str());
        return(simulationConfigCode);
      }
    }


    FunctionExit("Configure");

    return(0);
  }


  int application::ConfigureGeometries()
  {
    FunctionEntry("ConfigGeom");

    std::ostringstream localMessages;

    std::vector<std::string> geometryNames(appConfig.GetValueVector<std::string>("Geometry:Names"));

    int numGeometries = geometryNames.size();
    if(numGeometries == 0){
      ErrOut("Warning: No geometry configuration specified.\n");
      return(1);
    }
    geometryInfo.geometryNames = geometryNames;
    geometryInfo.gridNames.resize(numGeometries);
    geometryInfo.gridInfos.resize(numGeometries);
    for(int iGeo = 0;iGeo < numGeometries;iGeo++){
      std::string geometryName(geometryNames[iGeo]);
      std::vector<std::string> gridNames(appConfig.GetValueVector<std::string>(geometryName+":GridNames"));
      geometryInfo.gridNames[iGeo] = gridNames;
      int numGeomGrids = gridNames.size();
      geometryInfo.gridInfos[iGeo].resize(numGeomGrids);
      std::vector<simulation::geometry::gridinfo> &gridInfos(geometryInfo.gridInfos[iGeo]);
      for(int iGeomGrid = 0;iGeomGrid < numGeomGrids;iGeomGrid++){
        const std::string gridName(geometryName+":"+gridNames[iGeomGrid]);
        plascom2::util::ConfigureGridInfo(appConfig,gridName,gridInfos[iGeomGrid],
                                          localMessages);
      }
    }
    FunctionExit("ConfigGeom");

    return(0);
  };

  int application::ConfigureDomains()
  {

    FunctionEntry("ConfigDomains");

    std::ostringstream localMessages;

    std::vector<std::string> domainNames(appConfig.GetValueVector<std::string>("Domain:Names"));

    int numDomains = domainNames.size();

    // Should never do this because of check in caller
    if(numDomains == 0){
      numDomains = 1;
      domainNames.resize(numDomains,"Domain");
    }

    if(numDomains > 1){
      ErrOut("Error: Only single domain simulations are currently supported.\n");
      FunctionExit("ConfigDomains");
      return(1);
    }

    domainInfo.numDomains = numDomains;
    domainInfo.domainNames = domainNames;

    localMessages << "Number of domains: " << numDomains << std::endl
                  << "Domain Names: ";
    std::vector<std::string>::iterator dnIt = domainNames.begin();
    while(dnIt != domainNames.end())
      localMessages << *dnIt++ << " ";
    localMessages << std::endl;
    DebugOut(localMessages.str());
    pcpp::io::RenewStream(localMessages);

    // Not sure splitting should be done here, no harm when only one domain
    std::ostringstream messageStream;
    DebugOut("Splitting domains.\n");
    if(util::SplitDomains(appConfig,globalCommunicator,domainInfo,appDomains,messageStream)){
      ErrOut("Domain splitting failed:\n");
      ErrOut(messageStream.str());
      FunctionExit("ConfigDomains");
      return(1);
    } else {
      StdOut(messageStream.str(),2);
      pcpp::io::RenewStream(messageStream);
    }
    DebugOut("Domains split.\n");


    DebugOut("Configuring domains.\n");
    domainAdvancers.resize(numDomains,NULL);
    domainInitializers.resize(numDomains);
    for(int iDomain = 0;iDomain < numDomains;iDomain++){
      std::string &domainName(domainInfo.domainNames[iDomain]);
      localMessages << "Configuring domain: " << domainName << std::endl;
      DebugOut(localMessages.str());
      pcpp::io::RenewStream(localMessages);
      int globalDomainIndex = domainInfo.globalDomainIndex[iDomain];
      fixtures::ConfigurationType domainConfig(pcpp::util::ExtractConfigParams(appConfig,domainName));
      if(appDomains[iDomain].ConfigureDomain(domainConfig,messageStream)){
        messageStream << "Domain(" << iDomain+1 << ") \"" << domainName << "\""
                      << " failed to configure." << std::endl;
        ErrOut(messageStream.str());
        FunctionExit("ConfigDomains");
        return(1);
      }      
      // Set the Geometry Info 
      appDomains[iDomain].SetGeometryInfo(geometryInfo);
    }

    DebugOut("Domains configured.\n");

    FunctionExit("ConfigDomains");
    return(0);
  }

  std::string VirtualNodeReport(pcpp::VirtualNodeType &virtualNode)
  {
    std::ostringstream Ostr;
    Ostr << "VirtualNode:NumNodes       = " << virtualNode.numNodes       << std::endl
         << "VirtualNode:PPN            = " << virtualNode.procsPerNode   << std::endl
         << "VirtualNode:ioPerNode      = " << virtualNode.ioPerNode      << std::endl
         << "VirtualNode:performIO      = " << virtualNode.performIO      << std::endl
         << "VirtualNode:performCompute = " << virtualNode.performCompute << std::endl;
    if(virtualNode.ioRank > -1){
      Ostr << "VirtualNode:ioRank         = " << virtualNode.ioRank << "/" 
           << virtualNode.ioCommunicator.Size() << std::endl;
    }
    if(virtualNode.computeRank > -1) {
      Ostr << "VirtualNode:computeRank    = " << virtualNode.computeRank << "/"
           << virtualNode.computeCommunicator.Size() << std::endl;
    }
    return(Ostr.str());
  }


  int application::ConfigureSimulation()
  {
    FunctionEntry("ConfigSim");
    std::ostringstream messageStream; 

    // Configure the virtual node(s)
    int configureVirtualNodeCode = util::ConfigureVirtualNode(appConfig,globalCommunicator,
                                                              virtualNode,messageStream);

    // Report the virtual node confguration
    if(configureVirtualNodeCode){
      ErrOut("Configuration of virtual node failed.\n");
      return(configureVirtualNodeCode);
    } else {
      StdOut(messageStream.str());
      pcpp::io::RenewStream(messageStream);
    }
    StdOut(VirtualNodeReport(virtualNode),2);

    // Configure the geometries
    int configureGeometryCode = ConfigureGeometries();
    if(configureGeometryCode){
      ErrOut("Configuration of geometries failed.\n");
      return(configureGeometryCode);
    }

    // Configure the domains
    int configureDomainsCode = ConfigureDomains();
    if(configureDomainsCode){
      ErrOut("Configuration of domains failed.\n");
      return(configureDomainsCode);
    }
    
    numStepsMax       = appConfig.GetValue<int>(ConfigKey("PlasCom2","NumSteps"));
    if(numStepsMax < 0)
      numStepsMax = appConfig.GetValue<int>(ConfigKey("PlasCom2","NumStepsMax"));
    if(numStepsMax < 0)
      numStepsMax = 0;
    
    numStepsIO     = appConfig.GetValue<int>(ConfigKey("PlasCom2","NumStepsIO"));
    numStepsStatus = appConfig.GetValue<int>(ConfigKey("PlasCom2","NumStepsStatus"));
    

    // Make a list of all the grid names for the local domains 
    // domains might get split again according to grids.
    std::list<std::string> allLocalGridNames;
    int numLocalDomains = domainInfo.numLocalDomains;
    for(int iDomain = 0;iDomain < numLocalDomains;iDomain++){
      int globalDomainIndex = domainInfo.globalDomainIndex[iDomain];
      std::vector<std::string> domainGridNames(appDomains[iDomain].GridNames());
      std::vector<std::string>::iterator dnIt = domainGridNames.begin();
      while(dnIt != domainGridNames.end())
        allLocalGridNames.push_back(*dnIt++);
    }
    allLocalGridNames.sort();
    allLocalGridNames.unique();
    int numAllLocalGrids = allLocalGridNames.size();

    // ----- FIXME - right now hard-code to single domain/single grid -----
    if(numLocalDomains > 0){

      int globalDomainIndex = domainInfo.globalDomainIndex[0]; // hardcode to one (first) domain
      DomainBaseType &simDomain(appDomains[0]);
      DomainBaseType::gridvector &domainGrids(simDomain.Grids());
      DomainBaseType::statevector &gridStates(simDomain.States());
      DomainBaseType::statevector &gridParams(simDomain.Params());

      const std::vector<std::string> &domainGridNames(simDomain.GridNames());
      int numGlobalGrids = domainGridNames.size();

      pcpp::CommunicatorType &domainComm(simDomain.Communicator());
      int numProcDomain = domainComm.Size();
      int myDomainRank  = domainComm.Rank();

      // Partition the grids among processors
      simDomain.PartitionDomain(0);

      std::vector<int> &globalGridIndices(simDomain.LocalGridIndices());
      int numLocalGrids = simDomain.NumberOfLocalGrids();
      //      std::cout << "Number of local grids: " << numLocalGrids << std::endl;
      for(int iLocalGrid = 0;iLocalGrid < numLocalGrids;iLocalGrid++){
       
        int iGrid = globalGridIndices[iLocalGrid];
        //         std::cout << "Processing grid(" << iLocalGrid << "," << iGrid << ")" 
        //                   << std::endl;
        const std::string &gridName(domainGridNames[iGrid]);

        domainGrids[iGrid] = new GridType;
        gridStates[iGrid]  = new StateType;
        gridParams[iGrid]  = new StateType;

        GridType &domainGrid(*domainGrids[iGrid]);
        StateType &gridState(*gridStates[iGrid]);
        StateType &paramState(*gridParams[iGrid]);

        pcpp::CommunicatorType &domainGridComm(simDomain.GridCommunicator(iGrid));
        
        bool generateGrid = false;
        bool readGrid = false;
        std::string gridFileName;

        int gridConfigCode = plascom2::util::ConfigureGrid<GridType>(appConfig,gridName,domainGrid,messageStream); 

        if(gridConfigCode){
          messageStream << "Grid(" << iGrid << "," << iLocalGrid << "," << gridName
                        << ") failed to configure." << std::endl;
          ErrOut(messageStream.str());
          return(gridConfigCode);
        } else {
          StdOut("++++++++++++\n",2);
          StdOut("Grid Configuration Messages:\n",2);
          StdOut(messageStream.str(),2);
          StdOut("++++++++++++\n",2);
          pcpp::io::RenewStream(messageStream);
          if(appConfig.IsSet(ConfigKey(gridName,"GenerateGrid")))
            generateGrid = appConfig.GetFlagValue(ConfigKey(gridName,"GenerateGrid"));
          if(appConfig.IsSet(ConfigKey(gridName,"ReadGrid")))
            readGrid = appConfig.GetFlagValue(ConfigKey(gridName,"ReadGrid"));
          if(!generateGrid){
            if(appConfig.IsSet(ConfigKey(gridName,"GridFile"))){
              gridFileName = appConfig.GetValue(ConfigKey(gridName,"GridFile"));
            } else if(Restart()) {
              gridFileName = appConfig.GetValue("RestartFileName");
            } else {
              StdOut("WARNING: Setting up without grid data. Consequences unknown.\n");
            }
          }
          if(readGrid && gridFileName.empty()){
            ErrOut("ReadGrid set but no grid file specified.\n");
            return(1);
          }
        }
        
        // @todo Need to examine the following boundary handling block for multiple grids
        {
          //          std::cout << "Configuring grid domain boundaries for GridID: " << iGrid+1 << std::endl;
          int domainBoundaryConfigCode = simDomain.ConfigureGridDomainBoundaries(appConfig,iGrid,messageStream);
          
          if(domainBoundaryConfigCode){
            messageStream << "Domain boundaries failed to configure." << std::endl;
            ErrOut(messageStream.str());
            return(domainBoundaryConfigCode);
          } else {
            StdOut(messageStream.str(),2);
            pcpp::io::RenewStream(messageStream);
          }
        }

        int gridType = domainGrid.Type();
        messageStream << "Grid (" << gridName << ") has gridType(" 
                      << gridType << ")" << std::endl;
        if(true) {
          
          typedef GridType::HaloType HaloType;
          
          
          pcpp::CommunicatorType &gridComm(domainGrid.Communicator());      
          std::vector<size_t>    &gridSizes(domainGrid.GridSizes());
          std::vector<double>    &physicalExtent(domainGrid.PhysicalExtent());
          std::vector<bool>      &isPeriodic(domainGrid.PeriodicDirs());
          int numDim = gridSizes.size();

          if(isPeriodic.empty())
            isPeriodic.resize(numDim,false);

          
          size_t numGlobalNodes = 1;
          size_t numGlobalCells = 1;

          for(int iDim = 0;iDim < numDim;iDim++){
            numGlobalNodes *= gridSizes[iDim];
            numGlobalCells *= (gridSizes[iDim]-1);
          }
          
          if((gridType < simulation::grid::RECTILINEAR || generateGrid) && !readGrid){
            
            std::vector<double> dX(numDim,0);
            
            if(physicalExtent.empty()){
              physicalExtent.resize(numDim*2,0);
              for(int iDim = 0;iDim < numDim;iDim++)
                physicalExtent[2*iDim+1] = 1.0;
            }
            
            for(int iDim = 0;iDim < numDim;iDim++){
              dX[iDim] = (physicalExtent[2*iDim+1]-physicalExtent[2*iDim])/static_cast<double>(gridSizes[iDim]-1);
            }
            
            domainGrid.SetGridSpacings(dX);
          }
          
          int spatialOrder  = simDomain.SpatialOrder();

          if(!plascom2::operators::sbp::IsValidOrder(spatialOrder)){
            messageStream << "ERROR! Invalid spatial order specified." << std::endl;
            ErrOut(messageStream.str());
            return(1);
          }
          
          int interiorOrder  = (spatialOrder - 1)*2;
          DomainBaseType::OperatorType &inOperator(simDomain.Operator());
          plascom2::operators::sbp::Initialize(inOperator,interiorOrder);
          
          //        int boundaryDepth = inOperator.boundaryDepth;
          int boundaryDepth = std::max(3,inOperator.boundaryDepth+1);
          //          int boundaryDepth = std::max(3,interiorOrder+1);
          
          std::vector<int> haloDepths(2*numDim,boundaryDepth);
          
          //          pcpp::ParallelTopologyInfoType cartInfo;
          if(true){

            // determine which domain boundaries have bc's in this domain
            // Grab the grid-specific boundaries for this grid
            // Check the direction of each one
            std::vector<DomainBaseType::BoundaryType> &gridBoundaries(simDomain.GridBoundaries(iGrid));

            messageStream << "Processing " << gridBoundaries.size() << " boundaries." 
                          << std::endl;
            StdOut(messageStream.str(),2);
            pcpp::io::RenewStream(messageStream);
            
            std::vector<DomainBaseType::BoundaryType>::iterator boundIt = gridBoundaries.begin();
            
            int processingBoundary = 1;
            while(boundIt != gridBoundaries.end()){

              //              std::cout << "Boundary " << processingBoundary++ << std::endl;
              DomainBaseType::BoundaryType &gridBoundary(*boundIt++);
              DomainBaseType::BCType &boundaryCondition(gridBoundary.BC());
              DomainBaseType::GridType &boundaryGrid(gridBoundary.Grid());

              const std::string &boundaryName(gridBoundary.Name());
              std::vector<size_t> &boundaryGridSizes(boundaryGrid.GridSizes());
              assert(!boundaryGridSizes.empty());
              //              std::cout << "Boundary name: " << boundaryName << std::endl;
              std::vector<simulation::grid::subregion> &gridSubRegions(boundaryGrid.SubRegions());
              int numSubRegions = gridSubRegions.size();
              //              std::cout << "Number of subregions: " << numSubRegions << std::endl;
              int gridSubRegionIndex = gridBoundary.GridRegionID();
              if(gridSubRegionIndex >= numSubRegions){
                messageStream << "Boundary subregion out of range" << std::endl;
                return(1);
              }
              //              std::cout << "Subregion index: " << gridSubRegionIndex << std::endl;
              simulation::grid::subregion &gridRegion(gridSubRegions[gridSubRegionIndex]);
              int direction = std::abs(gridRegion.normalDirection)-1;
              pcpp::IndexIntervalType &regionInterval(gridRegion.regionInterval);
              assert(!regionInterval.empty());
              
              if(gridRegion.normalDirection != 0){
                //                std::cout << "Normal direction: " << gridRegion.normalDirection << std::endl;
                if((gridRegion.regionInterval[direction].first  == 0) || 
                   (gridRegion.regionInterval[direction].second == (boundaryGridSizes[direction] - 1))){
                  isPeriodic[direction] = false;
                } else {
                  //                  std::cout << "Lies off domain boundary." << std::endl;                  
                }
                messageStream << "Domain [" << boundaryName << "] boundary in direction "
                              << (gridRegion.normalDirection > 0 ? "-" : "+") << "X" 
                              << direction+1 << " has BC [" << boundaryCondition.BCName() 
                              << "]." << std::endl;
              } else {
                messageStream << "Domain [" << boundaryName << "] is a subregion of type ["
                              << boundaryCondition.BCName() << "]." << std::endl;
              }

            }

            messageStream << "Setting up domain with PERIODIC directions (";
            pcpp::io::DumpContents(messageStream,isPeriodic,",");
            messageStream << ")" << std::endl;
            StdOut(messageStream.str(),2);
            pcpp::io::RenewStream(messageStream);
          }
        
          pcpp::ParallelTopologyInfoType pbsCartInfo;
          pbsCartInfo.numDimensions = numDim;
          pbsCartInfo.isPeriodic.resize(numDim,0);
          for(int iDim = 0;iDim < numDim;iDim++)
            if(isPeriodic[iDim]) pbsCartInfo.isPeriodic[iDim] = 1; 
        
          // Check for grid-guidance on partitioning 
          // (any such guidance came from the grid config earlier)
          std::vector<int> &decompDirs(domainGrid.DecompDirs());
          if(!decompDirs.empty()){
            pbsCartInfo.cartDecompDirections = decompDirs;
          }
          
          pcpp::comm::SetupCartesianTopology(domainGridComm,pbsCartInfo,gridComm,messageStream);
        
          // Grab the local coordinates and global dimension of the Cartesian topo
          std::vector<int> &cartCoords(gridComm.CartCoordinates());
          std::vector<int> &cartDims(gridComm.CartDimensions());
          
          // Generate a report of the Cartesian setup and put it the messageStream
          pcpp::report::CartesianSetup(messageStream,pbsCartInfo);
        
          pcpp::IndexIntervalType &partitionInterval(domainGrid.PartitionInterval());
          pcpp::IndexIntervalType globalInterval;
          globalInterval.InitSimple(gridSizes);
        
          // === Partition the grid ===
          if(pcpp::util::PartitionCartesianInterval(globalInterval,cartDims,cartCoords,
                                                    partitionInterval,messageStream)){
            messageStream << "ConfigureDomain:SetupSimulationFixtures:Error: Partitioning failed." << std::endl;
            ErrOut(messageStream.str());
            return(1);
          }
        
          std::vector<size_t> extendGrid(2*numDim,0);
          //      std::vector<bool> isPeriodic(numDim,true);
          size_t numPointsPart = partitionInterval.NNodes();
        
          std::vector<bool> haveNeighbors(2*numDim,true);
          for(int iDim = 0;iDim < numDim; iDim++){
            if(isPeriodic[iDim]){
              extendGrid[2*iDim]   = haloDepths[2*iDim];
              extendGrid[2*iDim+1] = haloDepths[2*iDim+1];
            } else {
              if(partitionInterval[iDim].first == 0) {
                haveNeighbors[2*iDim] = false;
              } else { 
                extendGrid[2*iDim] = haloDepths[2*iDim];
              }
              if(partitionInterval[iDim].second == (gridSizes[iDim]-1)){
                haveNeighbors[2*iDim + 1] = false;
              } else {
                extendGrid[2*iDim + 1] = haloDepths[2*iDim+1];
              }
            }
          }
        
          messageStream << "ConfigureDomain: Configuring grid extensions (";
          pcpp::io::DumpContents(messageStream,extendGrid,",");
          messageStream << ")" << std::endl;
          StdOut(messageStream.str(),2);
          pcpp::io::RenewStream(messageStream);

          domainGrid.SetDimensionExtensions(extendGrid);

          if(domainGrid.Finalize()){ // generates coordinates if true is passed
            messageStream << "TestFixtures:SetupSimulationFixtures:Error: Grid finalization failed." << std::endl;
            return(1);
          }
          
          messageStream << simulation::report::Grid(domainGrid) << std::endl;
          StdOut(messageStream.str(),2);
          pcpp::io::RenewStream(messageStream);

        
          const std::vector<size_t> &bufferSizes(domainGrid.BufferSizes());
          //        const pcpp::IndexIntervalType &partitionInterval(domainGrid.PartitionInterval());
          const pcpp::IndexIntervalType &partitionBufferInterval(domainGrid.PartitionBufferInterval());  
          size_t numPointsBuffer = domainGrid.BufferSize();
      
          // Set up the HALO data structure for this grid/geometry 
          HaloType &inHalo(domainGrid.Halo());
          inHalo.SetGridInterval(globalInterval);
          inHalo.SetPartitionInterval(partitionInterval);

          // This forces the halo to have neighbors even on 
          // partitions with domain boundaries.
          inHalo.SetNeighbors(haveNeighbors);
          inHalo.SetPeriodicDirs(pbsCartInfo.isPeriodic);

          std::vector<pcpp::IndexIntervalType> remoteHaloExtents(inHalo.CreateRemoteHaloExtents(extendGrid));
          std::vector<pcpp::IndexIntervalType> localHaloExtents(inHalo.CreateLocalHaloExtents(extendGrid));
  
          inHalo.SetLocalBufferSizes(bufferSizes);
          inHalo.SetLocalPartitionExtent(partitionBufferInterval);
  
          inHalo.SetRemoteHaloExtents(remoteHaloExtents);
          inHalo.SetLocalHaloExtents(localHaloExtents);
  
          const std::vector<pcpp::IndexIntervalType> &remoteHaloBufferExtents(inHalo.RemoteHaloBufferExtents());
          const std::vector<pcpp::IndexIntervalType> &localHaloBufferExtents(inHalo.LocalHaloBufferExtents());
  
          inHalo.CreateSimpleSendIndices();
          inHalo.CreateSimpleRecvIndices();
        }

        if(true) {

          int numScalar = 0;
          bool withFlux = false;
          if(simDomain.CreateData(iGrid,messageStream)){
            messageStream << "Domain data creation failed." << std::endl;
            ErrOut(messageStream.str());
            pcpp::io::RenewStream(messageStream);
            return(1);
          } else {
            messageStream << "Domain data created." << std::endl;
            StdOut(messageStream.str(),2);
            pcpp::io::RenewStream(messageStream);
          }
        }
      
        int refLengthHandle = paramState.GetFieldIndex("refLength");
        if(refLengthHandle >= 0){
          double *refLengthPtr = paramState.GetFieldBuffer<double>("refLength");
          double refLengthValue = *refLengthPtr;
          if(refLengthValue > 0){
            std::vector<double> lengthScale(domainGrid.Dimension(),1.0/refLengthValue);
            domainGrid.SetScaleFactor(lengthScale);
          }
        }
      
        if(gridType < simulation::grid::RECTILINEAR && !readGrid){
          generateGrid = true;
        }
 
        bool coordinatesReady = false;
        int generatorReturnCode = 0;
        if(generateGrid && !readGrid){
          const std::string generatorMode(appConfig.GetValue(gridName+":GeneratorMode"));
          if(gridType < simulation::grid::RECTILINEAR || generatorMode == "Uniform" ||
             generatorMode == "Uniform-rectangular"){
            generatorReturnCode = domainGrid.GenerateCoordinates(messageStream);
          } else {
            generatorReturnCode = GenerateGrid(domainGrid,gridName,appConfig,messageStream);
          }
          if(generatorReturnCode == 0){
            coordinatesReady = true;
          } else {
            ErrOut("PC2Configuration::Generate grid coordinates failed with messages:\n");
            ErrOut(messageStream.str());
            return(1);
          }      
        } else if(!gridFileName.empty()){
          if(pcpp::io::hdf5::FetchGridData(domainGrid,gridName,gridFileName,messageStream)){
            ErrOut(std::string("Fetch grid data failed:\n")+messageStream.str()+std::string("\n"));
            return(1);
          } else {
            StdOut(std::string("Read grid data for (")+gridName+
                   std::string(") from ")+gridFileName+std::string(".\n"),2);
            StdOut(messageStream.str(),2);
            pcpp::io::RenewStream(messageStream);
            coordinatesReady = true;
          }
        } else {
          StdOut("WARNING: No grid coordinate data. Unsure what's next.");
        }
        domainGrid.SetType(gridType);
        
        if(coordinatesReady)
          domainGrid.ExchangeCoordinates();
        
      } // End of loop over local grids
      
      // Report the domain configuration
      //    int numDomains = appDomains.size();
      std::ostringstream domainReportStream;
      domainReportStream << "====== Domain Report =======" << std::endl;
      
      for(int iDomain = 0;iDomain < domainInfo.numLocalDomains;iDomain++){
        std::string domainName(domainInfo.domainNames[iDomain]);
        int globalDomainIndex2 = domainInfo.globalDomainIndex[iDomain];
        domainReportStream << "------" << domainName << "------" << std::endl;
        domainReportStream << "Domain ID (local,global) = (" << iDomain+1 << ","
                           << globalDomainIndex2+1 << ")" << std::endl;
        domainReportStream << appDomains[iDomain].Report() << std::endl;
      }
      
      StdOut(domainReportStream.str(),2);

#ifdef USE_OVERKIT
      int error = InitializeOverkit();
      if (error) {
        ErrOut("Overkit failed to initialize.\n");
        return(error);
      }
      ConfigureOverkit();
#endif
      
    }

    FunctionExit("ConfigSim");
    return(0);    
  }

#ifdef USE_OVERKIT
  int application::InitializeOverkit() {

    ovk_log_level LogLevel = OVK_LOG_ERRORS;
    if (VerbLevel() >= 1) LogLevel = ovk_log_level(LogLevel | OVK_LOG_WARNINGS);
    if (VerbLevel() >= 2) LogLevel = ovk_log_level(LogLevel | OVK_LOG_STATUS);

    ovk_context_params *contextParams;
    ovkCreateContextParams(&contextParams);
    ovkSetContextParamComm(contextParams, globalCommunicator.GetCommunicator());
    ovkSetContextParamLogLevel(contextParams, LogLevel);
    ovkSetContextParamErrorHandlerType(contextParams, OVK_ERROR_HANDLER_RETURN);

    int error = ovkCreateContext(&overkitData.context, contextParams);
    if (error != OVK_ERROR_NONE) {
      return(1);
    }

    ovkDestroyContextParams(&contextParams);

    return(0);

  }

  int application::FinalizeOverkit() {

    ovkDestroyContext(&overkitData.context);

    return(0);

  }

  int application::ConfigureOverkit() {

    overkitData.domains.resize(domainInfo.numLocalDomains);

    for (int iDomain = 0; iDomain < domainInfo.numLocalDomains; iDomain++) {

      DomainBaseType &simDomain(appDomains[iDomain]);
      DomainBaseType::gridvector &domainGrids(simDomain.Grids());
      const std::vector<std::string> &domainGridNames(simDomain.GridNames());

      // Overkit currently assumes all grids within a domain have the same dimension (now that
      // I think about it, I'm not sure this is a good assumption...)
      int numDim = domainGrids[0]->GridSizes().size();

      ovk_domain_params *domainParams;
      ovkCreateDomainParams(overkitData.context, &domainParams);
      ovkSetDomainParamName(domainParams, simDomain.Name().c_str());
      ovkSetDomainParamDimension(domainParams, numDim);
      ovkSetDomainParamComm(domainParams, simDomain.Communicator().GetCommunicator());

      ovkCreateDomain(overkitData.context, &overkitData.domains[iDomain], domainParams);

      ovkDestroyDomainParams(overkitData.context, &domainParams);

      DomainBaseType::gridvector::iterator gridIt = domainGrids.begin();

      for (int iGrid = 0; iGrid < simDomain.NumberOfGrids(); iGrid++) {

        int gridID = iGrid+1;

        // Assume single grid for now
//         if ((*gridIt)-><ID>() == gridID) {
        if (true) {

          GridType &domainGrid(**gridIt);

          pcpp::CommunicatorType  &gridComm(domainGrid.Communicator());
          pcpp::IndexIntervalType &partitionInterval(domainGrid.PartitionInterval());
          std::vector<bool>       &isPeriodic(domainGrid.PeriodicDirs());
          std::vector<int>        &cartCoords(gridComm.CartCoordinates());
          std::vector<int>        &cartDims(gridComm.CartDimensions());

          // Local index range
          int localBegin[3] = {0}, localEnd[3] = {1};
          for (int d = 0; d < numDim; ++d) {
            localBegin[d] = partitionInterval[d].first;
            localEnd[d] = partitionInterval[d].second + 1;
          }

          // Periodicity
          bool periodic[3] = {false};
          for (int d = 0; d < numDim; ++d) {
            periodic[d] = isPeriodic[d];
          }

          // Temporary hack to get the global grid size
          int numPointsGlobal[3] = {1};
          MPI_Allreduce(localEnd, numPointsGlobal, numDim, MPI_INT, MPI_MAX,
            gridComm.GetCommunicator());

          // Neighboring ranks
          int neighborOffsetLower[3] = {0}, neighborOffsetUpper[3] = {0};
          int numNeighbors = 1;
          for (int d = 0; d < numDim; ++d) {
            neighborOffsetLower[d] = cartCoords[d] > 0 ? -1 : 0;
            neighborOffsetUpper[d] = cartCoords[d] < cartDims[d]-1 ? 1 : 0;
            numNeighbors *= neighborOffsetUpper[d]-neighborOffsetLower[d]+1;
          }
          --numNeighbors; // Don't count self as a neighbor
          std::vector<int> neighborRanks(numNeighbors);
          int iNeighbor = 0;
          for (int k = neighborOffsetLower[2]; k <= neighborOffsetUpper[2]; ++k) {
            for (int j = neighborOffsetLower[1]; j <= neighborOffsetUpper[1]; ++j) {
              for (int i = neighborOffsetLower[0]; i <= neighborOffsetUpper[0]; ++i) {
                if (i != 0 || j != 0 || k != 0) {
                  int neighborOffset[3] = {i, j, k};
                  int neighborCoords[3] = {0};
                  for (int d = 0; d < numDim; ++d) {
                    neighborCoords[d] = cartCoords[d] + neighborOffset[d];
                  }
                  MPI_Cart_rank(gridComm.GetCommunicator(), neighborCoords,
                    &neighborRanks[iNeighbor]);
                  ++iNeighbor;
                }
              }
            }
          }

          ovk_grid_params *gridParams;
          ovkCreateGridParams(overkitData.domains[iDomain], &gridParams);
          ovkSetGridParamID(gridParams, gridID);
          ovkSetGridParamName(gridParams, domainGridNames[iGrid].c_str());
          ovkSetGridParamComm(gridParams, gridComm.GetCommunicator());
          ovkSetGridParamSize(gridParams, numPointsGlobal);
          ovkSetGridParamPeriodic(gridParams, periodic);
//           ovkSetGridParamPeriodicStorage(gridParams, <PeriodicStorage>);
//           ovkSetGridParamPeriodicLength(gridParams, <PeriodicLength>);
//           ovkSetGridParamGridType(gridParams, <GridType>);
          ovkSetGridParamLocalBegin(gridParams, localBegin);
          ovkSetGridParamLocalEnd(gridParams, localEnd);
          ovkSetGridParamNeighborRanks(gridParams, numNeighbors, &neighborRanks[0]);

          ovkCreateGridLocal(overkitData.domains[iDomain], NULL, gridParams);

          ovkDestroyGridParams(overkitData.domains[iDomain], &gridParams);

          ++gridIt;

        } else {

          ovkCreateGridRemote(overkitData.domains[iDomain], NULL);

        }

      }

    }

    return(0);
  }
#endif
  
}