RunPlasCom2.C

Go to the documentation of this file.

#include <iomanip>

#include "PlasCom2.H"
#include "EulerRHS.H"
#include "RK4Advancer.H"
#include "PCPPCommUtil.H"
#include "PCPPReport.H"
#include "PCPPIntervalUtils.H"
#include "Initialization.H"

void ReportTimings(std::vector<std::vector<double> > &inTimers,const char **timerNames,
                   unsigned int nTimes,std::ostream &outStream);

#ifdef USE_PAPI
void ReportCounters(long long **inCounters,const char **timerNames,
                    std::vector<std::string> &papiCounterNames,
                    unsigned int nTimes,std::ostream &outStream);
#endif

int plascom2::application::RunPlasCom2()
{
  
  //  typedef plascom2::operators::sbp::base  operator_t;
  typedef DomainBaseType::OperatorType    operator_t;
  //  typedef DomainBaseType::RHSType         rhs_t;

  typedef euler::rhs<DomainBaseType::GridType,DomainBaseType::StateType,operator_t> rhs_t;
  typedef simulation::grid::halo          halo_t;
  typedef pcpp::ParallelGlobalType        global_t;
  typedef pcpp::CommunicatorType          comm_t;
  typedef pcpp::IndexIntervalType         interval_t;

  //  typedef euler::rhs<grid_t,state_t,operator_t> rhs_t;

  std::ostringstream messageStream;

  int myRank       = globalCommunicator.Rank();
  numProcessors    = globalCommunicator.NProc();
  numThreads       = 1;

  //  operator_t sbpOperator;
  
  FunctionEntry("SimSetup");

  std::string &domainName(domainInfo.domainNames[0]);
  bool restartStep = Restart();  
  DomainBaseType &simDomain(appDomains[0]);

  DomainBaseType::gridvector  &simGrids(simDomain.Grids());
  DomainBaseType::statevector &gridStates(simDomain.States());
  DomainBaseType::statevector &gridParams(simDomain.Params());
  DomainBaseType::statevector &gridTargets(simDomain.Targets());

  operator_t &sbpOperator(simDomain.Operator());

  int numLocalGrids = simDomain.NumberOfLocalGrids();
  int spatialOrder  = simDomain.SpatialOrder();
  int interiorOrder  = (spatialOrder - 1)*2;

  const std::vector<std::string> &domainGridNames(simDomain.GridNames());

  if(numLocalGrids <= 0){
    ErrOut("Error: no local grids.\n");
    return(1);
  }
  std::vector<std::vector<int> > numThreadsGridDim(numLocalGrids);
  comm_t &domainCommunicator(simDomain.Communicator());
  int startStep = 0;
  double startTime = 0;

  for(int iLocalGrid = 0;iLocalGrid < numLocalGrids;iLocalGrid++){

    //    std::vector<int> &localGridIndices(simDomain.LocalGridIndices());

    int iGrid = simDomain.LocalGridIndex(iLocalGrid);

    grid_t     &simGrid(*simGrids[iGrid]);
    state_t    &simState(*gridStates[iGrid]);
    state_t    &paramState(*gridParams[iGrid]);
    state_t    &targetState(*gridTargets[iGrid]);
    
    const std::string &gridName(domainGridNames[iGrid]);
    
    StdOut(messageStream.str());
    pcpp::io::RenewStream(messageStream);
  
    std::vector<size_t> &gridSizes(simGrid.GridSizes());

    pcpp::IndexIntervalType &partitionInterval(simGrid.PartitionInterval());
    pcpp::IndexIntervalType &partitionBufferInterval(simGrid.PartitionBufferInterval());
    halo_t     &simHalo(simGrid.Halo());
    
    if(!Restart()){

      FunctionEntry("InitSoln");
      
      if(euler::util::InitializeQuiescentState(simGrid,simState,paramState,0,&messageStream)){
        ErrOut("Could not initialize flow. Messages:\n");
        ErrOut(messageStream.str());
        pcpp::io::RenewStream(messageStream);
        domainCommunicator.SetExit(1);
      } else {
        StdOut(messageStream.str());
        pcpp::io::RenewStream(messageStream);
      }
      
      if(appConfig.IsSet(ConfigKey(simDomain.Name(),"Init:CaseName"))){
        
        std::string initName(appConfig.GetValue(ConfigKey(simDomain.Name(),"Init:CaseName")));
        std::string initKey(ConfigKey(ConfigKey(simDomain.Name(),"Init"),initName));
        std::vector<double> initParams(appConfig.GetValueVector<double>(ConfigKey(initKey,"Params")));
        std::vector<int> initFlags(appConfig.GetValueVector<int>(ConfigKey(initKey,"Flags")));
        
      
        messageStream << "Initializing flow on Domain (" << simDomain.Name() 
                      << ") with initialization case (" << initName << ")." << std::endl
                      << "Parameters: (";
        pcpp::io::DumpContents(messageStream,initParams,",");
        messageStream << ")" << std::endl
                      << "Flags: (";
        pcpp::io::DumpContents(messageStream,initFlags,",");
        messageStream << ")" << std::endl;
        StdOut(messageStream.str(),2);
        pcpp::io::RenewStream(messageStream);
        
        if(InitializeSolution(simGrid,simState,paramState,simDomain.DomainBoundary(iGrid),
                              initName,initParams,initFlags,0,&messageStream)){
          ErrOut("Solution initialization failed. Messages:\n");
          ErrOut(messageStream.str());
          pcpp::io::RenewStream(messageStream);
          domainCommunicator.SetExit(1);
        } else {
          StdOut(messageStream.str());
          pcpp::io::RenewStream(messageStream);
        }
        
      }
      FunctionExit("InitSoln");
      
      if(domainCommunicator.Check())
        return(1);
      
    } else {
      FunctionEntry("RestartSoln");
      StdOut("Processing restart.\n");
      // Read the grid, state *and target state* from restart file
      if(!appConfig.IsSet("RestartFileName")){
        ErrOut("Restart file name not set. Aborting.\n");
        FunctionExit("RestartSoln");
        return(1);
      }
      std::string restartFileName(appConfig.GetValue("RestartFileName"));
      std::ostringstream messageOut;
      messageOut << "Restart File Name = " << restartFileName << std::endl;
      StdOut(messageOut.str());
      pcpp::io::RenewStream(messageOut);
      pcpp::io::simfileinfo restartFileInfo;
      if(myRank == 0){
        if(pcpp::io::hdf5::FileInfo(restartFileName,restartFileInfo,messageOut)){
          ErrOut(std::string("ERROR: FileInfo failed for ")+
                 restartFileName+std::string("\nMessages:\n")+
                 messageOut.str()+std::string("\n"));
          domainCommunicator.SetExit(1);
        }
      }
      if(domainCommunicator.Check())
        return(1);
      domainCommunicator.StreamBroadCast(restartFileInfo);

      if(restartFileInfo.simStep > 0)
        startStep = restartFileInfo.simStep;
      if(restartFileInfo.simTime > 0)
        startTime = restartFileInfo.simTime;
      
      std::ostringstream domainPathStream;
      if(restartFileInfo.formatBits.test(pcpp::io::ISLEGACY)){
        domainPathStream <<  "Group" << (iGrid < 9  ? "00" :
                                         iGrid < 99 ? "0"  :
                                         "") << iGrid+1;
        restartStep = false;
      } else {
        domainPathStream << domainName << "/grid" << iGrid+1;
      }
      std::string dataPath = domainPathStream.str();
      if(pcpp::io::hdf5::FetchStateData(restartFileName,dataPath,simGrid,simState,messageOut)){
        std::string fetchMessages(messageOut.str());
        pcpp::io::RenewStream(messageOut);
        messageOut << "Failed to get state data from " << restartFileName 
                   << ". Aborting with messages:" << std::endl
                   << fetchMessages << std::endl;
        ErrOut(messageOut.str());
        FunctionExit("RestartSoln");
        return(1);
      }
      
      FunctionExit("RestartSoln");
    }

    simState.InitializeFieldHandles();
    paramState.InitializeFieldHandles();
    targetState.InitializeFieldHandles();

    // Set up the target if needed
    std::string targetFileName;
    if(appConfig.IsSet("TargetFileName")){
      targetFileName = appConfig.GetValue("TargetFileName");
    } 
    if(!targetFileName.empty()){
      state_t readTarget(simState);
      std::ostringstream messageOut;
      pcpp::io::simfileinfo targetFileInfo;
      if(myRank == 0){
        if(pcpp::io::hdf5::FileInfo(targetFileName,targetFileInfo,messageOut)){
          ErrOut(std::string("ERROR: FileInfo failed for ")+
                 targetFileName+std::string("\nMessages:\n")+
                 messageOut.str()+std::string("\n"));
          domainCommunicator.SetExit(1);
        }
      }
      if(domainCommunicator.Check())
        return(1);
      domainCommunicator.StreamBroadCast(targetFileInfo);
      std::ostringstream domainPathStream;
      if(targetFileInfo.formatBits.test(pcpp::io::ISLEGACY)){
        domainPathStream <<  "Group" << (iGrid < 9  ? "00" :
                                         iGrid < 99 ? "0"  :
                                         "") << iGrid+1;
      } else {
        domainPathStream << domainName << "/grid" << iGrid+1;
      }
      std::string dataPath = domainPathStream.str();
      if(pcpp::io::hdf5::FetchStateData(targetFileName,dataPath,simGrid,readTarget,messageOut)){
        std::string fetchMessages(messageOut.str());
        pcpp::io::RenewStream(messageOut);
        messageOut << "Failed to get target data from " << targetFileName 
                   << ". Aborting with messages:" << std::endl
                   << fetchMessages << std::endl;
        ErrOut(messageOut.str());
        FunctionExit("RestartSoln");
        return(1);
      }
      targetState.CopyStateData(readTarget);
    } else {
      targetState.CopyStateData(simState);
    }
    
    // -- Set up the RHS --
    rhs_t *nsRHS = new rhs_t;  
    nsRHS->SetGlobal(*this);
    
    FunctionEntry("InitRHS");
    
    if(nsRHS->Initialize(simGrid,simState,paramState,targetState,sbpOperator)){
      return(1);
    }    

    //     FunctionEntry("SetDomainBCs");
    //     //  nsRHS.SetDomainBoundaries(simDomain.DomainBoundary(0));
    //     nsRHS.SetDomainBCs(simDomain.DomainBoundary(iGrid),simDomain.BCs());
    //     FunctionExit("SetDomainBCs");
    
    simDomain.SetRHS(iGrid,*nsRHS);

    FunctionExit("InitRHS");
  }
  
  FunctionEntry("SetDomainBCs");
  //  nsRHS.SetDomainBoundaries(simDomain.DomainBoundary(0));
  simDomain.SetDomainBCs();
  FunctionExit("SetDomainBCs");

  rk4advancer<domain_t> simAdvancer;

  FunctionEntry("InitAdvancer");
  simAdvancer.SetCommunication(false);
  simAdvancer.InitializeAdvancer(simDomain,*this,messageStream);
  FunctionExit("InitAdvancer");
  
  StdOut(messageStream.str(),2);
  pcpp::io::RenewStream(messageStream);
  
  
#ifdef USE_OMP
  numThreads = omp_get_max_threads();
#endif

#pragma omp parallel
  {
    int  myThreadId   = 0;
    bool masterThread = true;

#ifdef USE_OMP
    myThreadId   = omp_get_thread_num();
    masterThread = (myThreadId == 0);
#endif

// #pragma omp master 
//     {
//     //       SyncIO();
//     }
    
#pragma omp barrier
#pragma omp master
    {
      StdOut("Partitioning grids among threads.\n",3);
    }
    
    simDomain.InitThreads();
    
#pragma omp barrier
#pragma omp master
    {
      StdOut("Partitioning grid among threads done.\n",3);
    }
    
    
    if(errorState)
      ErrOut("Grid threading setup failed.\n");
    
#pragma omp barrier
#pragma omp master
    {
      StdOut("Threads set up, syncing advancer intervals.\n",3);
      FunctionEntry("SyncIntervals");
    }
    
#pragma omp barrier

    simAdvancer.SyncIntervals();

#pragma omp barrier
    
#pragma omp master
    {
      StdOut("Intervals sync'd, syncing states.\n",3);
      FunctionExit("SyncIntervals");
      FunctionEntry("SyncStates");
    }

#pragma omp barrier

    simAdvancer.SyncStates();
    
#pragma omp barrier
#pragma omp master
    {
      StdOut("States sync'd.\n",3);
      FunctionExit("SyncStates");
      FunctionExit("SimSetup");
    }

    simDomain.SetTime(startTime);

#pragma omp barrier
#pragma omp master
    {
      messageStream << "Running PlasCom2 using PBS grids."         << std::endl
                    << "DomainName:   " << domainName              << std::endl
                    << "SBP Order:    " << interiorOrder           << std::endl
                    << "Current Step: " << startStep               << std::endl
                    << "Current Time: " << startTime               << std::endl
                    << "NumStepsMax:  " << numStepsMax             << std::endl
                    << "NumStepsIO :  " << numStepsIO              << std::endl
                    << "NumStepsStat: " << numStepsStatus          << std::endl                    
                    << "Number of Local Grids:  " << numLocalGrids << std::endl;
      
      for(int iLocalGrid = 0;iLocalGrid < numLocalGrids;iLocalGrid++){

        int iGrid = simDomain.LocalGridIndex(iLocalGrid); 
        grid_t     &simGrid(*simGrids[iGrid]);
        const std::string &gridName(domainGridNames[iGrid]);
        const std::vector<size_t> &gridSizes(simGrid.GridSizes());
        pcpp::IndexIntervalType &partitionInterval(simGrid.PartitionInterval());
        pcpp::IndexIntervalType &partitionBufferInterval(simGrid.PartitionBufferInterval());
 
        messageStream << "GridName:     " << gridName             << std::endl
                      << "GridSizes:    (";
        pcpp::io::DumpContents(messageStream,gridSizes,",");
        messageStream << ")"                                   << std::endl
                      << "Partition: ";
        partitionInterval.PrettyPrint(messageStream);
        messageStream << std::endl
                      << "Buffer Sizes: (";
        pcpp::io::DumpContents(messageStream,simGrid.BufferSizes(),",");
        messageStream << ")" << std::endl
                      << "PBI:";
        partitionBufferInterval.PrettyPrint(messageStream);
        messageStream << std::endl
                      << "Number of Threads: (";
        pcpp::io::DumpContents(messageStream,numThreadsGridDim[iGrid],",");
        messageStream << ")" << std::endl;
        StdOut(messageStream.str(),2);
        pcpp::io::RenewStream(messageStream);
      }
      
    }
#pragma omp barrier
#pragma omp master
    {
      FunctionEntry("Stepping");
      StdOut("Beginning stepping.\n",3);
    }
    
#pragma omp barrier
    
    // ========= Time Stepping ===========
    //   errorState = 0;
    //   int iStep = 0;
    for(int iiStep = startStep;iiStep < numStepsMax;iiStep++){
      
#pragma omp master
      { 
        iStep = iiStep;
        simDomain.SetStep(iStep);
      }
      
#pragma omp barrier
      
      if(numStepsStatus > 0){
        if(!(iiStep%numStepsStatus)){
#pragma omp master
          {
            messageStream << " Step = "  << iiStep 
                          << " Time = "  << simDomain.Time() 
                          << std::endl;
            StdOut(messageStream.str(),1);
            pcpp::io::RenewStream(messageStream);
          }
        }
      }
      
      if(numStepsIO > 0){
        if(!(iiStep%numStepsIO) && !restartStep){
          
#pragma omp barrier
#pragma omp master
          {
            messageStream << "Computing DV before IO" << std::endl;
            StdOut(messageStream.str(),1);
            pcpp::io::RenewStream(messageStream);
          }
          
          simDomain.ComputeDV(myThreadId);
          
#pragma omp barrier
          
          if(errorState){
            ErrOut("ComputeDV failed.\n");
#pragma omp master
            {
              FunctionExit("Stepping");
            }
            break;
          }
          
#pragma omp master
          {
            messageStream << "Performing IO at step " << iiStep << std::endl;
            StdOut(messageStream.str(),1);
            pcpp::io::RenewStream(messageStream);
            FunctionEntry("IO");
            errorState = WriteDomains();
            FunctionExit("IO");
          }
          
#pragma omp barrier
          
          if(errorState){
            ErrOut("IO failed.\n");
#pragma omp master
            {
              FunctionExit("Stepping");
            }
            break;
          }
          
        }
      }
      restartStep = false;
#pragma omp barrier
#pragma omp master
      {      
        FunctionEntry("AdvanceDomain");
      }
      
      
      // This must be called for each domain being simulated
      int advanceCode = simAdvancer.AdvanceDomain();
      if(advanceCode){
        std::ostringstream errMsg;
        errMsg << "Advancer returned error code (" << advanceCode
               << "):" << std::endl << messageStream.str() 
               << std::endl;
        ErrOut(errMsg.str());
        errorState = 1;
      }
      
#pragma omp barrier
#pragma omp master
      {      
        FunctionExit("AdvanceDomain");
      }
      
      if(errorState)
        break;

      // if (domain.Status()) {
      //    Check on domain health and stopping condition
      // }
      
#pragma omp master
      {
        if(iiStep == 1) // this starts advancer timers after 2 steps
          simAdvancer.SetTimers(true);
      }
 
#pragma omp barrier
#pragma omp master
      {
        iStep++;
      }
    }
    
#pragma omp barrier
  }  // End OMP parllel section
  
  // Important: reset grid/rhs thread environment
  simDomain.SetNumThreads(1);
  simDomain.SetStep(iStep);

  int myThreadId = 0;

  if(errorState)
    return(1);
  
  FunctionExit("Stepping");
  
  messageStream << "Ending at step(" << iStep 
                << ") time( " << simDomain.Time() 
                << ")" << std::endl;

  StdOut(messageStream.str());
  pcpp::io::RenewStream(messageStream);
  
  // Check domain status to see if IO should be done
  messageStream << "Computing DV before final IO" << std::endl;
  StdOut(messageStream.str(),1);
  pcpp::io::RenewStream(messageStream);

  simDomain.ComputeDV(myThreadId);

  if(errorState){
    ErrOut("EulerPBSPeriodic: ComputeDV failed.\n");
  }

  messageStream << "Writing final dump." << std::endl;
  StdOut(messageStream.str(),1);
  pcpp::io::RenewStream(messageStream);

  FunctionEntry("IO");
  errorState = WriteDomains();
  FunctionExit("IO");

  if(errorState)
    return(1);
  
  // Set up for correct profiling output
  int numGrids = simDomain.NumberOfGrids();
  this->SetNumGrids(numGrids);
  for(int iGrid = 0;iGrid < numGrids;iGrid++){
    int gridDimension = simGrids[iGrid]->Dimension();
    size_t numPointsGrid = 1;
    std::vector<size_t> &gridSizes(simGrids[iGrid]->GridSizes());
    for(int iDim = 0;iDim < gridDimension;iDim++)
      numPointsGrid *= gridSizes[iDim];
    this->SetNumPointsGrid(iGrid,numPointsGrid);
  }
  this->SetNumSteps(iStep - startStep);

  double *processTimes;
  unsigned int nTimes;
  std::vector<std::vector<double> > parallelTimes;
  std::vector<std::vector<long long> > parallelCounts;
  const char **timerNames;
  simAdvancer.GetTimers(&processTimes,&timerNames,&nTimes);
  pcpp::comm::ReduceTimers(processTimes,nTimes,parallelTimes,globalCommunicator);
  if(myRank == 0){
    std::ostringstream profileStream;
    ReportTimings(parallelTimes,timerNames,nTimes,profileStream);
    StdOut(profileStream.str(),2);
    //     std::string profileFileName(profileName+processorTagOut.str());
    //     std::ofstream profileFile;
    //     profileFile.open(profileFileName.c_str());
    //     profileFile << profileStream.str() << std::endl;
    //     profileFile.close();
  }

    // #ifdef USE_PAPI
    //     long long **advancerCounters;
    //     testAdvancer.GetCounters(&advancerCounters);
    //     pcpp::io::RenewStream(profileStream);
    //     ReportCounters(advancerCounters,timerNames,papiCounterNames,nTimes,profileStream);
    //     std::string counterProfileName(std::string("PAPICounters")+processorTagOut.str());
    //     if(myRank == 0){
    //       std::ofstream counterFile;
    //       counterFile.open(counterProfileName.c_str());
    //       counterFile << profileStream.str() << std::endl;
    //       counterFile.close();
    //     }
    // #endif
    //   }
   

  return(0);
}


  
  
// #ifdef USE_PAPI
//   PAPI_library_init(PAPI_VER_CURRENT);
//   if(PAPI_create_eventset(&papiEventSet) != PAPI_OK){
//     PAPI_perror("PAPI event creation failed");
//     exit(1);
//   }
//   numCounters = 0;
//   std::vector<std::string>::iterator papiNameIt = papiEventNames.begin();
//   while(papiNameIt != papiEventNames.end()){
//     std::string eventName(*papiNameIt++);
//     int eventId = 0;
//     if(PAPI_event_name_to_code(const_cast<char *>(eventName.c_str()),&eventId) != PAPI_OK){
//       std::ostringstream papiErrStream;
//       papiErrStream << "Event not found (" << eventName << ")\n";
//       myGlobal.ErrOut(papiErrStream.str());
//       continue;
//     }
//     if(PAPI_add_event(papiEventSet,eventId) != PAPI_OK){
//       std::ostringstream papiErrStream;
//       papiErrStream << "Event could not be added (" << eventName << ")\n";
//       PAPI_perror(const_cast<char *>(papiErrStream.str().c_str()));
//       continue;
//     } else {
//       papiCounterNames.push_back(eventName);
//       numCounters++;
//     }
//   }
//   if(PAPI_start(papiEventSet) != PAPI_OK){
//     PAPI_perror("Could not start PAPI");
//     exit(1);
//   }
//   std::vector<long long> counterValues(numCounters,0);
//   //  std::vector<long long> counterTicks(numCounters,0);
// #endif
  
//   // Remove the previous probe file
//   pcpp::io::Remove("euler3d.probe");
  
//   int numThreads = 1;
//   int errorState = 0;
  
// #ifdef USE_OMP
//   numThreads = omp_get_max_threads();
// #endif
//   processorTagOut << "_"
//                   << (numThreads < 10   ? "000" :
//                       numThreads < 100  ? "00"  :
//                       numThreads < 1000 ? "0"   : "")
//                   << numThreads;
//   //  profileName = profileName + processorTagOut.str();
  

// #ifdef USE_PAPI
//  std::ostringstream papiDriverReport;
//  PAPI_stop(papiEventSet,&counterValues[0]);
//  papiDriverReport << "Driver PAPI Counters: "                     << std::endl;
//  for(int iCounter = 0;iCounter < numCounters;iCounter++)
//    papiDriverReport << papiCounterNames[iCounter] << ":\t\t" << counterValues[iCounter] << std::endl;
//  myGlobal.StdOut(papiDriverReport.str());
// #endif

//   messageStream << "Final Time: " << t << std::endl;

//   //  if(numStepsIO > 0){
//   //    domainState.Report(messageStream);
//   //    myGlobal.StdOut(messageStream.str());
//   //    pcpp::io::RenewStream(messageStream);
//   //  }
  
//   if(profLevel > 0){
//     std::ostringstream profileStream;
//     profileStream << "#Number of grids:   " << 1          << std::endl
//                   << "#Number of procs:   " << nProc      << std::endl
//                   << "#Number of threads: " << numThreads << std::endl
//                 << "#Number of points:  " << numGlobalNodes[0] << " " 
//                   << numGlobalNodes[1] << " " << numGlobalNodes[2] 
//                   << std::endl
//                   << "#Number of steps:   " << numSteps   << std::endl;
//     myGlobal.Report(profileStream);
//     double *processTimes;
//     unsigned int nTimes;
//     std::vector<std::vector<double> > parallelTimes;
//     std::vector<std::vector<long long> > parallelCounts;
//     const char **timerNames;
//     testAdvancer.GetTimers(&processTimes,&timerNames,&nTimes);
//     pcpp::comm::ReduceTimers(processTimes,nTimes,parallelTimes,simComm);
//     if(myRank == 0){
//       ReportTimings(parallelTimes,timerNames,nTimes,profileStream);
//       std::string profileFileName(profileName+processorTagOut.str());
//       std::ofstream profileFile;
//       profileFile.open(profileFileName.c_str());
//       profileFile << profileStream.str() << std::endl;
//       profileFile.close();
//     }
// #ifdef USE_PAPI
//     long long **advancerCounters;
//     testAdvancer.GetCounters(&advancerCounters);
//     pcpp::io::RenewStream(profileStream);
//     ReportCounters(advancerCounters,timerNames,papiCounterNames,nTimes,profileStream);
//     std::string counterProfileName(std::string("PAPICounters")+processorTagOut.str());
//     if(myRank == 0){
//       std::ofstream counterFile;
//       counterFile.open(counterProfileName.c_str());
//       counterFile << profileStream.str() << std::endl;
//       counterFile.close();
//     }
// #endif
//   }
  
//   myGlobal.Finalize();
  
//   return(0);

void ReportTimings(std::vector<std::vector<double> > &inTimers,const char **timerNames,unsigned int nTimes,
                   std::ostream &outStream)
{
  outStream << "# -------------------------------------------------------------------" << std::endl
            << "# Routine\tMin\t\tMax\t\tMean" << std::endl
            << "# ___________________________________________________________________" << std::endl;
  outStream.setf(std::ios::left);
  for(int iTimer = 0;iTimer < nTimes;iTimer++){
    outStream << std::setw(12) << timerNames[iTimer] 
              << "\t" << std::setw(14) << inTimers[0][iTimer]
              << "\t" << std::setw(14) << inTimers[1][iTimer] 
              << "\t" << std::setw(14) << inTimers[2][iTimer] 
              << std::endl;
  }
}

#ifdef USE_PAPI
void ReportCounters(long long **inCounters,const char **timerNames,std::vector<std::string> &papiCounterNames,
                    unsigned int nTimes,std::ostream &outStream)
{
  outStream << "# -----------------------------------------------------------------" << std::endl
            << "# Routine";
  for(int iCount = 0;iCount < numCounters;iCount++)
    outStream << "\t" << papiCounterNames[iCount];
  outStream << std::endl
            << "# _________________________________________________________________" << std::endl;
  outStream.setf(std::ios::left);
  for(int iTimer = 0;iTimer < nTimes;iTimer++){
    outStream << std::setw(12) << timerNames[iTimer];
    for(int iCount = 0;iCount < numCounters;iCount++){
      outStream << "\t" << std::setw(14) << inCounters[iTimer][iCount];
    }
    outStream << std::endl;
  }
}
#endif