PC2Util.C

Go to the documentation of this file.

#include <set>

#include "PC2Util.H"
#include "NavierStokesBC.H"
#include "NavierStokesDictionary.H"


using namespace simulation::domain::boundary;

namespace plascom2 {
  namespace util {

    inline int SymmetricElementID(int numDim,int i,int j)
    {
      int elementId = i*j;
      if(elementId <= numDim)
        return(elementId);
      int testId = elementId/numDim;
      if(testId < 2)
        return(elementId);
      if(testId == 2)
        return(elementId - 1);
      return(elementId - numDim);
    }

    int PC2Compare(const std::string &redFileName,const std::string &blueFileName,
                   double errTolerance,std::ostream &outStream)
    {
      
      std::ostringstream          errStream;
      pcpp::io::simfileinfo       redFileInfo;
      pcpp::io::simfileinfo       blueFileInfo;
      fixtures::ConfigurationType pc2Config;
      pcpp::CommunicatorType*     commPtr = NULL;

#ifdef ENABLE_HDF5
      if(pcpp::io::hdf5::FileInfo(redFileName,redFileInfo,errStream)){
        std::cerr << "ERROR: FileInfo failed for " << redFileName 
                  << ".\nHDF5 messages:" 
                  << errStream.str() << std::endl;
        return(-1);
      }
      if(pcpp::io::hdf5::FileInfo(blueFileName,blueFileInfo,errStream)){
        std::cerr << "ERROR: FileInfo failed for " << blueFileName 
                  << ".\nHDF5 messages:" 
                  << errStream.str() << std::endl;
        return(-1);
      }
#else
      return(-1);
#endif
      
      std::ostringstream infoStream;
      infoStream << "--------- File1 (Red) -----------" << std::endl;
      pcpp::report::SimFileInfo(infoStream,redFileInfo);
      infoStream << "--------- File2 (Blue) -----------" << std::endl;
      pcpp::report::SimFileInfo(infoStream,blueFileInfo);
      outStream << infoStream.str() << std::endl;      

      std::ostringstream processingStream;
      if(!pcpp::io::Compatible(redFileInfo,blueFileInfo,processingStream)){
        processingStream << "Error: Files are not compatible or cannot be processed." << std::endl;
        std::cerr << processingStream.str() << std::endl;
        return(-1);
      }

      int numGrids = redFileInfo.numGrids;
      for(int iGrid = 0;iGrid < numGrids;iGrid++){
        
        plascom2::grid_t  redGrid;
        plascom2::state_t redFileState;
        if(plascom2::io::hdf5::ReadSingle(redFileInfo,iGrid,redGrid,redFileState,processingStream)){
          processingStream << "Error: Could not get state and grid from " << redFileInfo.fileName
                           << std::endl;
          std::cerr << processingStream.str() << std::endl;
          return(-1);
        }
      
        // Debugging 
        //   std::ostringstream stateReportStream;
        //   redFileState.Report(stateReportStream);
        //   outStream << "Red State Report: " << std::endl
        //             << stateReportStream.str() << std::endl;
        //   pcpp::io::RenewStream(stateReportStream);
        
        plascom2::grid_t  blueGrid;
        plascom2::state_t blueFileState;
        if(plascom2::io::hdf5::ReadSingle(blueFileInfo,iGrid,blueGrid,blueFileState,processingStream)){
          processingStream << "Error: Could not get state and grid from " << blueFileInfo.fileName
                           << std::endl;
          std::cerr << processingStream.str() << std::endl;
          return(-1);
        }
      
        bool anyTest = false;
        bool gridTestPassed = true;
        
        
        std::vector<double> gridDiffNormData(4,0.0);
        int numDim        = 1;
        size_t numPoints  = 1;
        size_t numPoints2 = 1;
        numDim = redFileInfo.gridNumDimensions[0];
        const std::vector<double> &coordinateData1(redGrid.CoordinateData());  
        const std::vector<double> &coordinateData2(blueGrid.CoordinateData());
        numPoints  = coordinateData1.size()/numDim;
        numPoints2 = coordinateData2.size()/numDim;
        //   outStream << "Number of Points: " << numPoints << std::endl;
        if(numPoints  != numPoints2){
          outStream << "Incompatible grid sizes: Red(" << numPoints
                    << ") != Blue(" << numPoints2 << ")"
                    << std::endl;
          return(-1);
        }
        
        if(redFileInfo.formatBits.test(pcpp::io::HASGRID) &&
           blueFileInfo.formatBits.test(pcpp::io::HASGRID)){
          
          // Debugging
          //   blueFileState.Report(stateReportStream);
          //   outStream << "Blue State Report: " << std::endl
          //             << stateReportStream.str() << std::endl;
          //   pcpp::io::RenewStream(stateReportStream);
          
          std::vector<double> &normData(gridDiffNormData);
          
          anyTest = true;
          
          std::vector<double> gridDiff(numPoints*numDim,0);
          std::vector<double>::iterator gridDiffIt       = gridDiff.begin();
          std::vector<double>::const_iterator redGridIt  = coordinateData1.begin();
          std::vector<double>::const_iterator blueGridIt = coordinateData2.begin();
          while(gridDiffIt != gridDiff.end())
            *gridDiffIt++ = *redGridIt++ - *blueGridIt++;
          
          size_t maxErrLocation = 0;
          outStream << "# (Red Grid  - Blue Grid):" << std::endl;
          outStream << "# Coordinate    1-Norm     2-Norm    inf-Norm    Sum" << std::endl;
          for(int iDim = 0;iDim < numDim;iDim++){
            pcpp::util::ErrorMetrics(numDim,numPoints,&gridDiff[iDim*numPoints],&normData[0],maxErrLocation);
            outStream << iDim+1 << "     " << normData[1] << "    " <<  normData[2]
                      << "    " << normData[0] << "    " << normData[3] << std::endl;
            if(normData[0] > errTolerance)
              gridTestPassed = false;
            
          } 
        } else {
          outStream << "One or more grids not present, grid test skipped." << std::endl;
        }
        
        
        bool stateTestPassed = true;
        
        if(redFileInfo.formatBits.test(pcpp::io::HASSTATE) &&
           blueFileInfo.formatBits.test(pcpp::io::HASSTATE)){
          
          anyTest = true;
          
          std::vector<std::string> stateFields;
          
          stateFields.push_back("rho");
          stateFields.push_back("rhoV");
          stateFields.push_back("rhoE");
          stateFields.push_back("scalarVars");
          
          plascom2::state_t &redState(redFileState);
          plascom2::state_t &blueState(blueFileState);
          
          std::vector<double> redNormData(4,0.0);
          std::vector<double> blueNormData(4,0.0);
          std::vector<double> diffNormData(4,0.0);
          
          if(redFileInfo.formatBits.test(pcpp::io::ISLEGACY)){
            simulation::state::ConvertLegacyState(numDim,numPoints,redState);
          }
          if(blueFileInfo.formatBits.test(pcpp::io::ISLEGACY)){
            simulation::state::ConvertLegacyState(numDim,numPoints,blueState);
          }
          
          int numFields1 = redState.SetStateFields(stateFields);
          int numFields2 = blueState.SetStateFields(stateFields);
          
          std::ostringstream stateReportStream;
          redState.Report(stateReportStream);
          outStream << "Red State Report: " << std::endl
                    << stateReportStream.str() << std::endl;
          pcpp::io::RenewStream(stateReportStream);
          
          blueState.Report(stateReportStream);
          outStream << "Blue State Report: " << std::endl
                    << stateReportStream.str() << std::endl;
          pcpp::io::RenewStream(stateReportStream);
          
          std::ostringstream stateValidationStream;
          if(euler::util::ValidateState(redState,blueState,stateValidationStream)){
            outStream << "Red state invalid: " << std::endl
                      << stateValidationStream.str() 
                      << std::endl;
            return(-1);
          }
          if(euler::util::ValidateState(blueState,redState,stateValidationStream)){
            outStream << "Blue state invalid: " << std::endl
                      << stateValidationStream.str() 
                      << std::endl;
            return(-1);
          }
          
          {
            
            double *rhoErr  = redState.template GetFieldBuffer<double>("rho");
            double *rhoVErr = redState.template GetFieldBuffer<double>("rhoV");
            double *rhoEErr = redState.template GetFieldBuffer<double>("rhoE");
            
            int numScalars = 0;
            int scalarHandle = redState.GetDataIndex("scalarVars");
            double *scalarVarsErr = NULL;
            if(scalarHandle >= 0){
              scalarVarsErr = redState.template GetFieldBuffer<double>("scalarVars");
              numScalars = redState.Meta()[scalarHandle].ncomp;
            }
            //double *scalarPtr    = scalarData.Data<double>();
            //pcpp::field::dataset::DataBufferType &scalarData(inState.Field(scalarHandle));
            
            size_t maxErrLocation = 0;
            
            outStream << "# Red State Metrics:" << std::endl;
            outStream << "# Field     1-Norm     2-Norm     inf-Norm     Sum" << std::endl;      
            pcpp::util::ErrorMetrics(numDim,numPoints,rhoErr,&redNormData[0],maxErrLocation);
            outStream << "rho         " << redNormData[1] << "    " <<  redNormData[2]
                      << "    " << redNormData[0] << "    " << redNormData[3] << std::endl; 
            for(int iDim = 0;iDim < numDim;iDim++){
              pcpp::util::ErrorMetrics(numDim,numPoints,&rhoVErr[iDim*numPoints],&redNormData[0],maxErrLocation);
              outStream << "rhoV-" << iDim+1 << "     " << redNormData[1] << "    " << redNormData[2]
                        << "    " << redNormData[0] << "    " << redNormData[3] << std::endl;
            } 
            pcpp::util::ErrorMetrics(numDim,numPoints,rhoEErr,&redNormData[0],maxErrLocation);
            outStream << "rhoE        " << redNormData[1] << "    " <<  redNormData[2]
                      << "    " << redNormData[0] << "    " << redNormData[3] << std::endl; 
            for(int iScalar = 0;iScalar < numScalars;iScalar++){
              pcpp::util::ErrorMetrics(numScalars,numPoints,&scalarVarsErr[iScalar*numPoints],
                                       &redNormData[0],maxErrLocation);
              outStream << "scalarVars-" << iScalar+1 << "     " << redNormData[1] 
                        << "    " << redNormData[2] << "    " << redNormData[0] 
                        << "    " << redNormData[3] << std::endl;
            } 
            
            rhoErr  = blueState.template GetFieldBuffer<double>("rho");
            rhoVErr = blueState.template GetFieldBuffer<double>("rhoV");
            rhoEErr = blueState.template GetFieldBuffer<double>("rhoE");
            numScalars = 0;
            scalarHandle = blueState.GetDataIndex("scalarVars");
            scalarVarsErr = NULL;
            if(scalarHandle >= 0){
              scalarVarsErr = blueState.template GetFieldBuffer<double>("scalarVars");
              numScalars = blueState.Meta()[scalarHandle].ncomp;
            }
            
            outStream << std::endl << std::endl 
                      << "# Blue State Metrics:" << std::endl;
            outStream << "# Field     1-Norm     2-Norm     inf-Norm     Sum" << std::endl;      
            pcpp::util::ErrorMetrics(numDim,numPoints,rhoErr,&blueNormData[0],maxErrLocation);
            outStream << "rho         " << blueNormData[1] << "    " <<  blueNormData[2]
                      << "    " << blueNormData[0] << "    " << blueNormData[3] << std::endl; 
            for(int iDim = 0;iDim < numDim;iDim++){
              pcpp::util::ErrorMetrics(numDim,numPoints,&rhoVErr[iDim*numPoints],
                                       &blueNormData[0],maxErrLocation);
              outStream << "rhoV-" << iDim+1 << "     " << blueNormData[1] << "    " << blueNormData[2]
                        << "    " << blueNormData[0] << "    " << blueNormData[3] << std::endl;
            } 
            
            pcpp::util::ErrorMetrics(numDim,numPoints,rhoEErr,&blueNormData[0],maxErrLocation);
            outStream << "rhoE        " << blueNormData[1] << "    " <<  blueNormData[2]
                      << "    " << blueNormData[0] << "    " << blueNormData[3] << std::endl; 
            for(int iScalar = 0;iScalar < numScalars;iScalar++){
              pcpp::util::ErrorMetrics(numScalars,numPoints,&scalarVarsErr[iScalar*numPoints],
                                       &blueNormData[0],maxErrLocation);
              outStream << "scalarVars-" << iScalar+1 << "     " << blueNormData[1] 
                        << "    " << blueNormData[2] << "    " << blueNormData[0] 
                        << "    " << blueNormData[3] << std::endl;
            } 
          }
          
          plascom2::state_t diffState(redState);
          diffState -= blueState;
          
          
          double *rhoErr  = diffState.template GetFieldBuffer<double>("rho");
          double *rhoVErr = diffState.template GetFieldBuffer<double>("rhoV");
          double *rhoEErr = diffState.template GetFieldBuffer<double>("rhoE");
          int numScalars = 0;
          int scalarHandle = diffState.GetDataIndex("scalarVars");
          double *scalarVarsErr = NULL;
          if(scalarHandle >= 0){
            scalarVarsErr = diffState.template GetFieldBuffer<double>("scalarVars");
            numScalars = diffState.Meta()[scalarHandle].ncomp;
          }
          
          size_t maxErrLocation = 0;
          
          outStream << std::endl;
          outStream << "# (Red - Blue) State Metrics (Absolute): " << std::endl;
          outStream << "# Field     1-Norm     2-Norm     inf-Norm     Sum" << std::endl;      
          pcpp::util::ErrorMetrics(numDim,numPoints,rhoErr,&diffNormData[0],maxErrLocation);
          outStream << "rho         " << diffNormData[1] << "    " <<  diffNormData[2]
                    << "    " << diffNormData[0] << "    " << diffNormData[3] << std::endl;
          if(std::abs(blueNormData[0]) > 0.0){
            if(std::abs(diffNormData[0])/std::abs(blueNormData[0]) > errTolerance ||
               std::isnan(diffNormData[0]) || std::isinf(diffNormData[0])) stateTestPassed = false;
          } else {
            if(std::abs(diffNormData[0]) > errTolerance ||
               std::isnan(diffNormData[0]) || std::isinf(diffNormData[0])) stateTestPassed = false;
          }
          
          for(int iDim = 0;iDim < numDim;iDim++){
            pcpp::util::ErrorMetrics(numDim,numPoints,&rhoVErr[iDim*numPoints],
                                     &diffNormData[0],maxErrLocation);
            outStream << "rhoV-" << iDim+1 << "     " << diffNormData[1] << "    " << diffNormData[2]
                      << "    " << diffNormData[0] << "    " << diffNormData[3] << std::endl;
            
            if(std::abs(blueNormData[0]) > 0.0){
              if(std::abs(diffNormData[0])/std::abs(blueNormData[0]) > errTolerance ||
                 std::isnan(diffNormData[0]) || std::isinf(diffNormData[0])) stateTestPassed = false;
            } else {
              if(std::abs(diffNormData[0]) > errTolerance ||
                 std::isnan(diffNormData[0]) || std::isinf(diffNormData[0])) stateTestPassed = false;
            }
          }
          
          
          pcpp::util::ErrorMetrics(numDim,numPoints,rhoEErr,&diffNormData[0],maxErrLocation);
          outStream << "rhoE        " << diffNormData[1] << "    " <<  diffNormData[2]
                    << "    " << diffNormData[0] << "    " << diffNormData[3] << std::endl; 
          
          if(std::abs(blueNormData[0]) > 0.0){
            if(std::abs(diffNormData[0])/std::abs(blueNormData[0]) > errTolerance ||
               std::isnan(diffNormData[0]) || std::isinf(diffNormData[0])) stateTestPassed = false;
          } else {
            if(std::abs(diffNormData[0]) > errTolerance ||
             std::isnan(diffNormData[0]) || std::isinf(diffNormData[0])) stateTestPassed = false;
          }
          
          for(int iScalar = 0;iScalar < numScalars;iScalar++){
            pcpp::util::ErrorMetrics(numScalars,numPoints,&scalarVarsErr[iScalar*numPoints],
                                     &diffNormData[0],maxErrLocation);
            outStream << "scalarVars-" << iScalar+1 << "     " << diffNormData[1] 
                      << "    " << diffNormData[2] << "    " << diffNormData[0] 
                      << "    " << diffNormData[3] << std::endl;
            
            if(std::abs(blueNormData[0]) > 0.0){
              if(std::abs(diffNormData[0])/std::abs(blueNormData[0]) > errTolerance ||
                 std::isnan(diffNormData[0]) || std::isinf(diffNormData[0])) stateTestPassed = false;
            } else {
              if(std::abs(diffNormData[0]) > errTolerance ||
                 std::isnan(diffNormData[0]) || std::isinf(diffNormData[0])) stateTestPassed = false;
            }
          }
          
        } else {
          outStream << "One or more states not present, state test skipped." << std::endl;
        }
      
        if(!stateTestPassed)
          outStream << "Two states were not equivalent within tolerance." << std::endl;
        if(!gridTestPassed)
          outStream << "Two grids were not equivalent within tolerance."  << std::endl;
        
        if(!stateTestPassed || !gridTestPassed || !anyTest)
          return(1);
      }
      
      return(0);
    };
    
    bool operator==(const legacybc &lhs,const legacybc &rhs)
    {
      if(lhs.gridID != rhs.gridID ||
         lhs.bcType != rhs.bcType ||
         lhs.bcDir  != rhs.bcDir)
        return(false);
      for(int iVal = 0;iVal < 6;iVal++)
        if(lhs.bcInterval[iVal] != rhs.bcInterval[iVal])
          return(false);
      return(true);
    };
    bool operator!=(const legacybc &lhs,const legacybc &rhs)
    {
      return(!(lhs==rhs));
    };
    std::ostream &operator<<(std::ostream &outStream,const legacybc &bcData)
    {
      outStream << bcData.gridID << " " << bcData.bcType << " " << bcData.bcDir
                << " " << bcData.bcInterval[0] << " " << bcData.bcInterval[1]
                << " " << bcData.bcInterval[2] << " " << bcData.bcInterval[3]
                << " " << bcData.bcInterval[4] << " " << bcData.bcInterval[5];
      return(outStream);
    };

    std::istream &operator>>(std::istream &inStream,legacybc &bcData)
    {
      inStream >> bcData.gridID >> bcData.bcType >> bcData.bcDir
               >> bcData.bcInterval[0] >> bcData.bcInterval[1]
               >> bcData.bcInterval[2] >> bcData.bcInterval[3]
               >> bcData.bcInterval[4] >> bcData.bcInterval[5];
      return(inStream);
    };
          
    int ReadLegacyBCDat(const std::string &bcFileName,std::vector<legacybc> &bcDat,
                        fixtures::CommunicatorType &inCommunicator,std::ostream &messageStream)
    {
      int myRank = inCommunicator.Rank();
      int returnStatus = 0;
      std::string bcString;
      if(myRank == 0){
        std::ifstream Inf;
        std::ostringstream Ostr;
        Inf.open(bcFileName.c_str());
        if(!Inf){
          messageStream << "Error: Could not open legacy BC data file ("
                        << bcFileName << "). Aborting read." << std::endl;
          inCommunicator.SetErr(1);
        } else {
          std::string bcLine;
          while(std::getline(Inf,bcLine)){
            std::string::size_type x = bcLine.find("#");
            if(x == 0)
              continue;
            if(x != std::string::npos){
              bcLine = bcLine.substr(0,x);
            }
            if(bcLine.empty())
              continue;
            std::istringstream Istr(bcLine);
            legacybc myBc;
            Istr >> myBc;
            Ostr << myBc << std::endl;
          }
          bcString = Ostr.str();
        }
      }
      if(inCommunicator.Check())
        returnStatus = 1;
      if(!returnStatus){
        inCommunicator.BroadCast(bcString);
        std::string bcLine;
        std::istringstream Istr(bcString);
        while(std::getline(Istr,bcLine)){
          legacybc myBC;
          std::istringstream bcInStr(bcLine);
          bcInStr >> myBC;
          bcDat.push_back(myBC);
        }
      }
      return(returnStatus);
    }

    // # shock capturing
    // # SHOCK_CAPTURE = TRUE : Updated Kawai & Lele (See Kawai, Shankar, Lele, JCP, 2010)
    // SHOCK_CAPTURE = FALSE
    // SHOCK_CAPTURE_HyperCmu = 0.002
    // SHOCK_CAPTURE_HyperCbeta = 1.75
    // SHOCK_CAPTURE_HyperCKappa = 0.01
    // SHOCK_CAPTURE_CLIP = .FALSE.
    // # sponge zone
    // SPONGE_AMP = 10.
    // SPONGE_POW = 2
    
    int ConvertLegacyConfiguration(fixtures::ConfigurationType &inConfig,
                                   fixtures::CommunicatorType &inCommunicator,
                                   std::ostream &messageStream)
    {
      int myRank = inCommunicator.Rank();
      
      
      fixtures::ConfigurationType::ParamType inParam;

      inConfig.SetParameter("PlasCom2:RunMode","PC2");
      inConfig.SetParameter("Domain:Names","PlasComCM");

      const std::string geomName("cmgeom");
      inConfig.SetParameter("Geometry:Names",geomName);

      double bcWallTemp = 298.0;
      double tempRef    = 298.0;
      if(inConfig.IsSet("TEMPERATURE_REFERENCE"))
        tempRef = inConfig.GetValue<double>("TEMPERATURE_REFERENCE");
      if(inConfig.IsSet("BCIC_WALL_TEMP"))
        bcWallTemp = inConfig.GetValue<double>("BCIC_WALL_TEMP");
      else
        bcWallTemp = tempRef;
      
      inConfig.ResetKey("NSTEPMAX","PlasCom2:NumSteps");
      inConfig.ResetKey("NRESTART","PlasCom2:NumStepsIO");
      inConfig.ResetKey("NOUTPUT","PlasCom2:NumStepsStatus");
      // inConfig.SetParameter("PlasCom2:NumSteps","100");
      // inConfig.SetParameter("PlasCom2:NumStepsStatus","1");
      // inConfig.SetParameter("PlasCom2:NumStepsIO","100");
      
      inConfig.ResetKey("SPACEORDER","PlasComCM:SpatialOrder");
      inConfig.ResetKey("LENGTH_REFERENCE","PlasComCM:Param:refLength");
      inConfig.ResetKey("TEMPERATURE_REFERENCE","PlasComCM:Param:refTemperature");
      inConfig.ResetKey("DENSITY_REFERENCE","PlasComCM:Param:refRho");
      inConfig.ResetKey("PRESSURE_REFERENCE","PlasComCM:Param:refPressure");
      inConfig.ResetKey("GAMMA_REFERENCE","PlasComCM:Param:gamma");
      inConfig.ResetKey("REYNOLDS_NUMBER","PlasComCM:Param:refRe");
      inConfig.ResetKey("PRANDTL_NUMBER","PlasComCM:Param:refPr");
      inConfig.ResetKey("TIMESTEP","PlasComCM:Param:inputDT");

      // artificial dissipation setup
      bool   artDiss     = false;
      if(inConfig.IsSet("SAT_ARTIFICIAL_DISSIPATION"))
        artDiss = inConfig.GetFlagValue("SAT_ARTIFICIAL_DISSIPATION");
      double adSigma = 0.0;
      double adDil   = 0.0;
      double adDilCutOff = 0.0;
      if(artDiss){
        adSigma     = 30.0;
        adDil       = 30.0;
        adDilCutOff = -1.5;
        if(inConfig.IsSet("AMOUNT_SAT_ARTIFICIAL_DISSIPATION"))
          adSigma = inConfig.GetValue<double>("AMOUNT_SAT_ARTIFICIAL_DISSIPATION");
        if(inConfig.IsSet("AMOUNT_SAT_ARTIFICIAL_DISSIPATION_FOR_DILATATION"))
          adDil = inConfig.GetValue<double>("AMOUNT_SAT_ARTIFICIAL_DISSIPATION_FOR_DILATATION");
        if(inConfig.IsSet("SAT_ARTIFICIAL_DISSIPATION_DIL_CUTOFF"))
          adDilCutOff = inConfig.GetValue<double>("SAT_ARTIFICIAL_DISSIPATION_DIL_CUTOFF");
      }
      inConfig.SetParameter("PlasComCM:Param:sigmaDissipation",adSigma);
      inConfig.SetParameter("PlasComCM:Param:sigmaDilatation",adDil);
      inConfig.SetParameter("PlasComCM:Param:dilatationCutoff",adDilCutOff);

      int numDim = inConfig.GetValue<int>("ND");
      int numStepsStatus = inConfig.GetValue<int>("PlasCom2:NumStepsIO");
      if(numStepsStatus >= 100)
        numStepsStatus /= 5;
      else
        numStepsStatus = 10;
      inConfig.SetParameter<int>("PlasCom2:NumStepsStatus",numStepsStatus);
      int numScalars = 0;
      if(inConfig.IsSet("NUMBER_OF_SCALARS"))
        numScalars = inConfig.GetValue<int>("NUMBER_OF_SCALARS");
      if(numScalars < 0) numScalars = 0;

      inConfig.SetParameter<double>("PlasComCM:Param:power",0.666);
      inConfig.SetParameter<double>("PlasComCM:Param:beta",1.0);
      inConfig.SetParameter<double>("PlasComCM:Param:bulkViscFac",0.6);
      inConfig.SetParameter<int>("PlasComCM:Param:Flags",1);
      inConfig.SetParameter<int>("PlasComCM:Param:nonDimensional",2);
      std::string bcFileName(inConfig.GetValue("BC_FNAME"));
      if(bcFileName.empty())
        bcFileName = "bc.dat";

      std::vector<legacybc> bcData;
      std::map<int,int> legacyBCToPC2;
      legacyBCToPC2[21] = bc::navierstokes::SAT_SLIP_ADIABATIC;
      legacyBCToPC2[22] = bc::navierstokes::SAT_NOSLIP_ISOTHERMAL;
      legacyBCToPC2[24] = bc::navierstokes::SAT_FARFIELD;
      legacyBCToPC2[99] = bc::navierstokes::SPONGE;
            
      messageStream << "Reading legacy boundary conditions from " << bcFileName
                    << "." << std::endl;

      if(ReadLegacyBCDat(bcFileName,bcData,inCommunicator,messageStream)){
        messageStream << "Read of legacy bc data failed." << std::endl;
        return(1);
      }
      int numBC = bcData.size();
      std::vector<std::string> gridNames;

      std::multimap<int,int> gridIdToLegacyBCIndex;
      std::multimap<int,std::string> bcToRegion;
      std::set<int> gridIds;
      for(int iBCLine = 0;iBCLine < numBC;iBCLine++){
        legacybc &legacyBC(bcData[iBCLine]);
        gridIdToLegacyBCIndex.insert(std::make_pair(legacyBC.gridID,iBCLine));
        gridIds.insert(legacyBC.gridID);
      }

      std::ostringstream gridNamesStream;
      std::ostringstream geomGridNamesStream;
      int numLegacyGrids = gridIds.size();
      messageStream << "Found configuration for " << numLegacyGrids 
                    << " legacy grids." << std::endl;
      for(int iGrid = 0;iGrid < numLegacyGrids;iGrid++){
        int gridId = iGrid+1;
        std::ostringstream nameStream;
        nameStream << "Group"
                   << (gridId >= 100 ? ""   :
                       gridId >= 10  ? "0"  :
                       "00") << gridId;
        const std::string gridName(nameStream.str());
        pcpp::io::RenewStream(nameStream);
        gridNames.push_back(gridName);
        gridNamesStream << gridName << " ";
        nameStream << geomName << ":" << gridName;
        const std::string geomGridName(nameStream.str());
        geomGridNamesStream << geomGridName << " ";
        pcpp::io::RenewStream(nameStream);

        // get the legacy BCs for this grid
        std::pair<std::multimap<int,int>::iterator,
          std::multimap<int,int>::iterator> gridBCs =
          gridIdToLegacyBCIndex.equal_range(gridId);
        int regionId = 1;

        // loop over bc lines for this grid
        while(gridBCs.first != gridBCs.second){
          int bcIndex = gridBCs.first->second;
          legacybc &legacyBC(bcData[bcIndex]);
          gridBCs.first++;

          messageStream << "legacyBC("<< bcIndex << "):" << std::endl
                        << "GridID: " << legacyBC.gridID << std::endl
                        << "BCType: " << legacyBC.bcType << std::endl
                        << "BCDir:  " << legacyBC.bcDir  << std::endl
                        << "BCInterval: (";
          for(int iDim = 0;iDim < numDim;iDim++){
            if(iDim > 0) messageStream << "x";
            messageStream << "[" << legacyBC.bcInterval[iDim*2] << ":"
                          << legacyBC.bcInterval[iDim*2+1] << "]";
          }
          


          std::ostringstream regionNameStream;
          regionNameStream << "Region" << regionId++;
          const std::string regionName(regionNameStream.str());
          pcpp::io::RenewStream(regionNameStream);
          nameStream << regionName << " ";
          regionNameStream << geomName << ":" << gridName
                           << ":" << regionName;
          std::string regionPath(regionNameStream.str());  
        
          int pc2BCID = legacyBCToPC2[legacyBC.bcType];
          bcToRegion.insert(std::make_pair(pc2BCID,regionPath));
          std::ostringstream regionConfigStream;
          regionConfigStream << legacyBC.bcDir << " ";
          for(int iDir = 0;iDir < 6;iDir++)
            regionConfigStream << legacyBC.bcInterval[iDir] << " ";
          inConfig.SetParameter(regionPath,regionConfigStream.str());

        } // loop over grid-specific bc lines
        inConfig.SetParameter(geomGridName+":RegionNames",nameStream.str());
        pcpp::io::RenewStream(nameStream);

      } // Loop over grids
      
      
      std::string geomGridParam(geomName+std::string(":GridNames"));
      inConfig.SetParameter(geomGridParam,gridNamesStream.str());
      inConfig.SetParameter("PlasComCM:GridNames",geomGridNamesStream.str());

      std::string gridFileName;
      pcpp::io::simfileinfo gridFileInfo;
      if(inConfig.IsSet("GRID_FNAME")){
        gridFileName = inConfig.GetValue("GRID_FNAME");
        if(myRank == 0 && ix::sys::FILEEXISTS(gridFileName)){
          std::ostringstream errStream;
          if(pcpp::io::hdf5::FileInfo(gridFileName,gridFileInfo,errStream)){
            messageStream << "ERROR: FileInfo failed for grid file (" 
                          << gridFileName << "). HDF5 messages:" << std::endl
                          << errStream.str() << std::endl;
            inCommunicator.SetErr(1);
          } 
        } 
        if(inCommunicator.Check())
          return(1);
        inCommunicator.StreamBroadCast(gridFileInfo); 
      }
      
      for(int iGrid = 0;iGrid < numLegacyGrids;iGrid++){
        const std::string &gridName(gridNames[iGrid]);
        std::string gridParam(geomName+":"+gridName);
        std::ostringstream paramNameStream;
        paramNameStream << gridParam << ":Dimension";
        inConfig.SetParameter(paramNameStream.str(),numDim);
        pcpp::io::RenewStream(paramNameStream);
        if(!gridFileName.empty()){
          paramNameStream << gridParam << ":GridFile";
          inConfig.SetParameter(paramNameStream.str(),gridFileName);
          pcpp::io::RenewStream(paramNameStream);
        }
        if(gridFileInfo.numGrids > 0){
          paramNameStream << gridParam << ":NumPoints";
          const std::vector<size_t> &gridSize(gridFileInfo.gridSizes[iGrid]);
          std::ostringstream gridSizeStream;
          for(int iDim = 0;iDim < numDim;iDim++)
            gridSizeStream << gridSize[iDim] << " ";
          inConfig.SetParameter(paramNameStream.str(),gridSizeStream.str());
          pcpp::io::RenewStream(paramNameStream);
        }
        paramNameStream << gridParam << ":GenerateGrid";
        inConfig.SetParameter(paramNameStream.str(),"NO");
        pcpp::io::RenewStream(paramNameStream);
      }


      // bc configuration setup
      double satSigma1   = 1.0;
      if(inConfig.IsSet("AMOUNT_SAT_SIGMAI1"))
        satSigma1 = inConfig.GetValue<double>("AMOUNT_SAT_SIGMAI1");
      double satSigma2   = 1.0;
      if(inConfig.IsSet("AMOUNT_SAT_SIGMAI2"))
        satSigma2 = inConfig.GetValue<double>("AMOUNT_SAT_SIGMAI2");

      double farfieldSigma1 =  .5;
      double farfieldSigma2 = 1.0;
      if(inConfig.IsSet("AMOUNT_SAT_SIGMAI1_FARFIELD"))
        farfieldSigma1 = inConfig.GetValue<double>("AMOUNT_SAT_SIGMAI1_FARFIELD");
      if(inConfig.IsSet("AMOUNT_SAT_SIGMAI2_FARFIELD"))
        farfieldSigma2 = inConfig.GetValue<double>("AMOUNT_SAT_SIGMAI2_FARFIELD");

      double spongeAmp   = 5.0;
      double spongePower = 2.0;
      if(inConfig.IsSet("SPONGE_AMP"))
        spongeAmp = inConfig.GetValue<double>("SPONGE_AMP");
      if(inConfig.IsSet("SPONGE_POW"))
        spongePower = inConfig.GetValue<double>("SPONGE_POW");

      std::ostringstream satSigmaStream;
      satSigmaStream << satSigma1 << " " << satSigma2;
      const std::string satSigma(satSigmaStream.str());

      std::ostringstream ffSigmaStream;
      ffSigmaStream << farfieldSigma1 << " " << farfieldSigma2
                    << " 0.0";
      const std::string ffSigma(ffSigmaStream.str());

      std::ostringstream noSlipStream;
      noSlipStream << satSigma1 << " " << satSigma2 << " "
                   << bcWallTemp << " 0.0";
      const std::string noSlipSetting(noSlipStream.str());
      
      std::ostringstream spongeStream;
      spongeStream << spongeAmp << " " << spongePower;
      const std::string spongeSetting(spongeStream.str());
      
      std::ostringstream bcNamesStream;
      std::multimap<int,std::string>::iterator bc2RegionsIt = bcToRegion.begin();
      while(bc2RegionsIt != bcToRegion.end()){
        int currentBC = bc2RegionsIt->first;
        int bcId = currentBC;
        std::ostringstream bcRegionStream;
        while(bcId == currentBC){
          bcRegionStream << bc2RegionsIt++->second << " ";
          bcId = bc2RegionsIt->first;
        }
        const std::string bcRegion(bcRegionStream.str());

        if(currentBC == bc::navierstokes::SAT_SLIP_ADIABATIC){
          bcNamesStream << "SlipWall ";
          inConfig.SetParameter("PlasComCM:SlipWall:BCType","SAT_SLIP_ADIABATIC");
          inConfig.SetParameter("PlasComCM:SlipWall:Param:Fields","Numbers Flags");
          inConfig.SetParameter("PlasComCM:SlipWall:Field:Numbers:Meta","2 s 8");
          inConfig.SetParameter("PlasComCM:SlipWall:Field:Flags:Meta","1 s 4");
          inConfig.SetParameter("PlasComCM:SlipWall:Field:Numbers",satSigma);
          inConfig.SetParameter("PlasComCM:SlipWall:Param:Flags","0");
          inConfig.SetParameter("PlasComCM:SlipWall:RegionNames",bcRegion);
        } else if (currentBC == bc::navierstokes::SAT_FARFIELD){
          bcNamesStream << "Farfield ";
          inConfig.SetParameter("PlasComCM:Farfield:BCType","SAT_FARFIELD");
          inConfig.SetParameter("PlasComCM:Farfield:Param:Fields","Numbers Flags");
          inConfig.SetParameter("PlasComCM:Farfield:Field:Numbers:Meta","3 s 8");
          inConfig.SetParameter("PlasComCM:Farfield:Field:Flags:Meta","1 s 4");
          inConfig.SetParameter("PlasComCM:Farfield:Param:Numbers",ffSigma);
          inConfig.SetParameter("PlasComCM:Farfield:Param:Flags","0");
          inConfig.SetParameter("PlasComCM:Farfield:RegionNames",bcRegion);
        } else if (currentBC == bc::navierstokes::SAT_NOSLIP_ISOTHERMAL){
          bcNamesStream << "NoSlip ";
          inConfig.SetParameter("PlasComCM:NoSlip:BCType","SAT_NOSLIP_ISOTHERMAL");
          inConfig.SetParameter("PlasComCM:NoSlip:Param:Fields","Numbers Flags");
          inConfig.SetParameter("PlasComCM:NoSlip:Field:Numbers:Meta","4 s 8");
          inConfig.SetParameter("PlasComCM:NoSlip:Field:Flags:Meta","1 s 4");
          inConfig.SetParameter("PlasComCM:NoSlip:Param:Numbers",noSlipSetting);
          inConfig.SetParameter("PlasComCM:NoSlip:Param:Flags","0");
          inConfig.SetParameter("PlasComCM:NoSlip:RegionNames",bcRegion);
        } else if(currentBC == bc::navierstokes::SPONGE){
          bcNamesStream << "Sponge ";
          inConfig.SetParameter("PlasComCM:Sponge:BCType","SPONGE");
          inConfig.SetParameter("PlasComCM:Sponge:Param:Fields","Numbers Flags");
          inConfig.SetParameter("PlasComCM:Sponge:Field:Numbers:Meta","2 s 8");
          inConfig.SetParameter("PlasComCM:Sponge:Field:Flags:Meta","1 s 4");
          inConfig.SetParameter("PlasComCM:Sponge:Param:Numbers",spongeSetting);
          inConfig.SetParameter("PlasComCM:Sponge:Param:Flags","0");
          inConfig.SetParameter("PlasComCM:Sponge:RegionNames",bcRegion);
        }
        pcpp::io::RenewStream(bcRegionStream);
      }
      inConfig.SetParameter("PlasComCM:BCNames",bcNamesStream.str());

      
      // This generates the required data dictionary to run a NavierStokes simulation
      bool writeDV = inConfig.GetFlagValue("WRITE_DV");
      bool writeMetrics = inConfig.GetFlagValue("WRITE_MET");

      
      navierstokes::CreateDictionaryConfiguration(numDim,3,numScalars,writeDV,true,true,writeMetrics,"PlasComCM",inConfig);
     
      
      return(0);
    }

    int ConfigureGridInfo(const fixtures::ConfigurationType &inConfig,
                          const std::string &gridName,
                          simulation::geometry::gridinfo &gridInfo,
                          std::ostream &messageStream)
    {
      messageStream << "Configuring Grid Inforomation for [" << gridName << "]" << std::endl;
      std::vector<size_t> gridSizes(inConfig.GetValueVector<size_t>(ConfigKey(gridName,"NumPoints")));
      int numDim = gridSizes.size();
      gridInfo.gridSizes = gridSizes;
      int gridType = simulation::grid::UNIRECT;
      std::string defaultGridTypeName("Uniform-rectangular");
      if(inConfig.IsSet(ConfigKey(gridName,"Type"))){
        std::string typeName(inConfig.GetValue(ConfigKey(gridName,"Type")));
        if(typeName == "Uniform") typeName = defaultGridTypeName;
        int gType = simulation::grid::ResolveTopoName(typeName);
        if(gType >= 0 && gType < simulation::grid::NUMTOPOTYPE)
          gridType = gType;
        else
          messageStream << "Warning: Invalid grid type specified (" << typeName 
                        << "). Defaulting to grid type (" << defaultGridTypeName
                        << ")." << std::endl;
      }
      std::string typeName(simulation::grid::ResolveTopoType(gridType));
      messageStream << "Grid (" << gridName << ") Type: " << typeName << std::endl;
      gridInfo.gridType = gridType;
      if(inConfig.IsSet(ConfigKey(gridName,"PhysicalExtent"))){
        std::vector<double> physicalExtent(inConfig.GetValueVector<double>(ConfigKey(gridName,"PhysicalExtent")));
        gridInfo.physicalExtent = physicalExtent;
      }
      if(inConfig.IsSet(ConfigKey(gridName,"DecompDirs"))){
        std::vector<int> decompDirs(inConfig.GetValueVector<int>(ConfigKey(gridName,"DecompDirs")));
        gridInfo.decompDirs = decompDirs;
      }
      if(inConfig.IsSet(ConfigKey(gridName,"DecompSizes"))){
        std::vector<int> decompSizes(inConfig.GetValueVector<int>(ConfigKey(gridName,"DecompSizes")));
        gridInfo.decompSizes = decompSizes;
      }
      if(inConfig.IsSet(ConfigKey(gridName,"ThreadDecompDirs"))){
        std::vector<int> threadDecompDirs(inConfig.GetValueVector<int>(ConfigKey(gridName,"ThreadDecompDirs")));
        gridInfo.threadDecompDirs = threadDecompDirs;
      }
      if(inConfig.IsSet(ConfigKey(gridName,"PeriodicDirs"))){
        std::vector<int> periodicDirs(inConfig.GetValueVector<int>(ConfigKey(gridName,"PeriodicDirs")));
        gridInfo.periodicDirs = periodicDirs;
      }
      if(inConfig.IsSet(ConfigKey(gridName,"PeriodicLengths"))){
        std::vector<double> periodicLengths(inConfig.GetValueVector<double>(ConfigKey(gridName,"PeriodicLengths")));
        gridInfo.periodicLengths = periodicLengths;
      }
      if(inConfig.IsSet(ConfigKey(gridName,"FileName"))){
        gridInfo.fileName = inConfig.GetValue(ConfigKey(gridName,"FileName"));
      }
      return(0);
    };




  }
}