RunMaxwellSolver.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"
#include "MaxwellSolver.H"
#include "MaxwellUtil.H"

#define EPS0            8.85418782e-12  // permittivity of free space [F/m]
#define MU0             1.25663706e-6   // permeability of free space [H/m]



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::RunMaxwellSolver()
{
  
  //  typedef plascom2::operators::sbp::base  operator_t;
  typedef DomainBaseType::OperatorType    operator_t;
  typedef simulation::grid::halo          halo_t;
  typedef pcpp::ParallelGlobalType        global_t;
  typedef pcpp::CommunicatorType          comm_t;
  typedef pcpp::IndexIntervalType         interval_t;

  typedef Maxwell::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]);

  DomainBaseType &simDomain(appDomains[0]);
  DomainBaseType::gridvector &simGrids(simDomain.Grids());
  DomainBaseType::statevector &gridStates(simDomain.States());
  DomainBaseType::statevector &gridParams(simDomain.Params());

  operator_t &sbpOperator(simDomain.Operator());

  //  std::vector<grid_t *> &simGrids(simDomain.Grids());
  grid_t     &simGrid(*simGrids[0]);
  state_t    &simState(*gridStates[0]);
  state_t    &paramState(*gridParams[0]);
  
  //  grid_t simGrid;
  //  state_t simState;
  //  state_t paramState;

  comm_t &domainCommunicator(simDomain.Communicator());

  const std::vector<std::string> &domainGridNames(simDomain.GridNames());
  const std::string &gridName(domainGridNames[0]);

  int spatialOrder  = simDomain.SpatialOrder();
  int interiorOrder = (spatialOrder - 1)*2;
  
  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)) {
//      ErrOut("Could not initialize flow.\n");
//      domainCommunicator.SetExit(1);
//    }

    /* ======== Initialize in entire buffer ======== */
    std::vector<double> EfieldAll(3, 0.0);
    std::vector<double> BfieldAll(3, 0.0);
    std::vector<double> JAll(3, 0.0);
    double epsilonAll = EPS0;
    double muAll      = MU0;

    if(Maxwell::util::InitializeMaxwellStateConstFields(simGrid, simState, &EfieldAll[0], &BfieldAll[0], &JAll[0], epsilonAll, muAll, true))
    {
      ErrOut("Failed to initialize in entire buffer.\n");
      domainCommunicator.SetExit(1);
    }

    /* ======== Initialize in partition only ======== */
//    std::vector<double> Efield(3,  1.0); Efield[1] = -2.0; Efield[2] =  3.0;
//    std::vector<double> Bfield(3, -1.0); Bfield[1] =  2.0; Bfield[2] = -3.0;
//    std::vector<double> J(3, 4.0);       J[1]      = -5.0; J[2]      =  6.0;
//    double epsilon = EPS0;
//    double mu      = MU0;
//
//    if(Maxwell::util::InitializeMaxwellStateConstFields(simGrid, simState, &Efield[0], &Bfield[0], &J[0], epsilon, mu, false))
//    {
//      ErrOut("Failed to initialize in partition.\n");
//      domainCommunicator.SetExit(1);
//    }

    /* ======== Initialize parameters ======== */
    double *CFLPtr  = paramState.template GetFieldBuffer<double>("CFL");
    double CFLValue = *CFLPtr;

    if(Maxwell::util::InitializeMaxwellParameters(simGrid, paramState, CFLValue))
    {
      ErrOut("Failed to initialize parameters.\n");
      domainCommunicator.SetExit(1);
    }


    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(0),
                            initName,initParams,initFlags,0)) {
        ErrOut("Solution initialization failed.\n");
        domainCommunicator.SetExit(1);
      }

    }
    FunctionExit("InitSoln");

    if(domainCommunicator.Check())
      return(1);

  }

  /* ======== Set up the RHS ======== */
  rhs_t MaxwellRHS;
  
  FunctionEntry("InitRHS");

  if(MaxwellRHS.Initialize(simGrid,simState,paramState,sbpOperator)) {
    return(1);
  }

  FunctionEntry("SetDomainBoundaries");
  MaxwellRHS.SetDomainBoundaries(simDomain.DomainBoundary(0));
  MaxwellRHS.SetDomainBCs(simDomain.BCs());
  FunctionExit("SetDomainBoundaries");
  simDomain.SetRHS(MaxwellRHS);

  FunctionExit("InitRHS");
  


  rk4advancer<domain_t> simAdvancer;

  FunctionEntry("InitAdvancer");
  simAdvancer.InitializeAdvancer(simDomain,*this,messageStream);
  FunctionExit("InitAdvancer");
  
  StdOut(messageStream.str(),2);
  pcpp::io::RenewStream(messageStream);
  
  
  // BEGIN OMP PARALLEL SECTION HERE
#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

    //     if(masterThread)
    //       SyncIO();
    
    //
    //
    
    
#pragma omp barrier

    if(masterThread) {
      StdOut("Partitioning grid among threads.\n",3);
    }
    
    // Insert thread partitioning here 
    // - initialize RHS and grid thread intervals
    std::vector<int> numThreadsDim;
    if(simGrid.SetupThreads(numThreadsDim))
      errorState = 1;
    MaxwellRHS.InitThreadIntervals();

#pragma omp barrier

    if(masterThread) {
      StdOut("Partitioning grid among threads done.\n",3);
    }


    if(errorState)
      ErrOut("Grid threading setup failed.\n");

#pragma omp barrier

    if(masterThread) {
      StdOut("Threads set up, syncing advancer intervals.\n",3);
    }

    if(masterThread)
      FunctionEntry("SyncIntervals");

#pragma omp barrier

    simAdvancer.SyncIntervals();

#pragma omp barrier
    
    if(masterThread) {
      StdOut("Intervals sync'd, syncing states.\n",3);
    }

    if(masterThread) {
      FunctionExit("SyncIntervals");
      FunctionEntry("SyncStates");
    }
#pragma omp barrier

    simAdvancer.SyncStates();
    
#pragma omp barrier

    if(masterThread) {
      StdOut("States sync'd.\n",3);
    }

    if(masterThread) {
      FunctionExit("SyncStates");
      FunctionExit("SimSetup");
    }

#pragma omp barrier
    
    if(masterThread) {
      messageStream << "Running PlasCom2 with PBS grid with: "  << std::endl
                    << "DomainName:   " << domainName           << std::endl
                    << "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,numThreadsDim,",");
      messageStream << ")" << std::endl << std::endl
                    << "SBP Order:    " << interiorOrder     << std::endl
                    << "Current Step: " << iStep             << std::endl
                    << "NumStepsMax:  " << numStepsMax       << std::endl
                    << "NumStepsIO :  " << numStepsIO        << std::endl
                    << "NumStepsStat: " << numStepsStatus    << std::endl;
      StdOut(messageStream.str(),2);
      pcpp::io::RenewStream(messageStream);
    }

    
#pragma omp barrier
    
    if(masterThread) {
      FunctionEntry("Stepping");
      StdOut("Beginning stepping.\n",3);
    }

#pragma omp barrier

    /* ======== Time Stepping ======== */
    //   errorState = 0;
    //   int iStep = 0;
    for(int iiStep = 0; iiStep < numStepsMax; iiStep++) {
      
      if(masterThread)
        iStep = iiStep;

#pragma omp barrier

      if(numStepsStatus > 0) {
        if(!(iiStep%numStepsStatus)) {
          if(masterThread) {
            messageStream << "Step = " << iiStep << std::endl;
            StdOut(messageStream.str(),1);
            pcpp::io::RenewStream(messageStream);
          }
        }
      }

      if(numStepsIO > 0) {
        if(!(iiStep%numStepsIO)) {

#pragma omp barrier

//          if(masterThread) {
//            messageStream << "Computing DV before IO" << std::endl;
//            StdOut(messageStream.str(),1);
//            pcpp::io::RenewStream(messageStream);
//          }
          
//          errorState = eulerRHS.ComputeDV(myThreadId);
          
//#pragma omp barrier
          
//          if(errorState) {
//            ErrOut("ComputeDV failed.\n");
//            if(masterThread) {
//              FunctionExit("Stepping");
//            }
//            break;
//          }

          if(masterThread) {
            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");
            if(masterThread) {
              FunctionExit("Stepping");
            }
            break;
          }
          
        }
      }
      
#pragma omp barrier
      
      if(masterThread) {
        FunctionEntry("AdvanceDomain");
      }
      
      
      // This must be called for each domain being simulated
      if(simAdvancer.AdvanceDomain()) {
        ErrOut("Advancer failed:\n");
        ErrOut(messageStream.str());
        errorState = 1;
      }
      
#pragma omp barrier
      
      if(masterThread) {
        FunctionExit("AdvanceDomain");
      }
      
      if(errorState)
        break;

#pragma omp barrier

      if(masterThread) {
        messageStream << "Computing DV" << std::endl;
        StdOut(messageStream.str(),1);
        pcpp::io::RenewStream(messageStream);
      }

      errorState = MaxwellRHS.ComputeDV(myThreadId);

#pragma omp barrier

      if(errorState) {
        ErrOut("ComputeDV failed.\n");
        if(masterThread) {
          FunctionExit("Stepping");
        }
        break;
      }


      // if (domain.Status()) {
      //    Check on domain health and stopping condition
      // }
      
      if(iiStep == 1) // this starts advancer timers after 2 steps
        simAdvancer.SetTimers(true);
      
 #pragma omp barrier

      if(masterThread)
        iStep++;

    }
    
#pragma omp barrier
  }  // End OMP parallel section
  
  // Important: reset grid/rhs thread environment
  simGrid.SetNumThreads(1);
  MaxwellRHS.InitThreadIntervals();

  int myThreadId = 0;

  if(errorState)
    return(1);
  
  FunctionExit("Stepping");
  
  messageStream << "Ending at step " << iStep << 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);

//  FunctionEntry("IO");
//  errorState = eulerRHS.ComputeDV(myThreadId);
//  FunctionExit("IO");

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

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

  if(errorState)
    return(1);
  
  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