RK4Advancer.H

Go to the documentation of this file.

#include "Advancer.H"
#include "PCPPTypes.H"
#include "StateOperators.H"
#include "State.H"
#include "OperatorKernels.H"
#include "RKKernels.H"
#include <limits>
//#include "RKUtil.H"
#ifdef USE_OMP
#include <omp.h>
#endif

//#define USE_PAPI 1
#ifdef USE_PAPI
#include <papi.h>
int papiEventSet = PAPI_NULL;
int numCounters;
#define MAXCOUNTERS 8

inline void AccumulateCounters(long long *tick,long long *tock,long long *clock)
{
  for(int iCounter = 0;iCounter < numCounters;iCounter++)
    clock[iCounter] += (tock[iCounter] - tick[iCounter]);
};

#endif

#ifndef _ADVANCER_TIMERS_
#define _ADVANCER_TIMERS_
static const char *timerNames[] = { "ZAXPY","RKSUM" };
#endif

namespace error {
  static const int INIT  = 1;
  static const int RHS   = 2;
  static const int IGX   = 4;
  static const int FINAL = 8;
}


//template<typename DomainType,typename RHSType>
template<typename DomainType>
class rk4advancer : public simulation::advancer::base<DomainType> {
  
public:

  enum {ZAXPYKERNEL=0,RKSUMKERNEL,NUMTIMERS};
  
  
  typedef simulation::advancer::base<DomainType> DomainAdvancerType;

  using DomainAdvancerType::messageStreamPtr;
  using DomainAdvancerType::domainPtr;
  using DomainAdvancerType::advancerConfig;
  
  //  using typename DomainType::GridType;

  typedef typename DomainType::GridType          GridType;
  typedef typename DomainType::StateType         StateType;
  typedef typename DomainType::HaloType          HaloType;
  typedef typename DomainType::OperatorType      OperatorType;
  typedef typename DomainType::RHSType           RHSType;

  bool errCheck;

  rk4advancer() : DomainAdvancerType(), // rkStage(0), 
                  numStages(4), rkH(0.0)
  {
    stageWeights.resize(4,0.5);
    stageWeights[2] = 1.0;
    stageWeights[3] = 1.0/6.0;
    updateWeights.resize(4,1.0/3.0);
    updateWeights[0] = 1.0/6.0;
    updateWeights[3] = 1.0/6.0;
    timersOn = false;
    communicateState = false;
  };

  rk4advancer(DomainType &inDomain) : DomainAdvancerType(inDomain), // rkStage(0),
                                      numStages(4), rkH(0.0) 
  {
    stageWeights.resize(4,0.5);
    stageWeights[2] = 1.0;
    stageWeights[3] = 1.0/6.0;
    updateWeights.resize(4,1.0/3.0);
    updateWeights[0] = 1.0/6.0;
    updateWeights[3] = 1.0/6.0;
    InitilializeAdvancer(inDomain,std::cout);
    timersOn = false;
    communicateState = false;
  };

  rk4advancer(fixtures::ConfigurationType &inConfig,
              DomainType &inDomain) : DomainAdvancerType(inDomain), // rkStage(0), 
                                      numStages(4), rkH(0.0) 
  {
    stageWeights.resize(4,0.5);
    stageWeights[2] = 1.0;
    stageWeights[3] = 1.0/6.0;
    updateWeights.resize(4,1.0/3.0);
    updateWeights[0] = 1.0/6.0;
    updateWeights[3] = 1.0/6.0;
    timersOn = false;
    communicateState = false;
    Configure(inConfig);
    InitilializeAdvancer(inDomain,std::cout);
  };

  void SetCommunication(bool onoff)
  { communicateState = onoff; };
  void SetTimers(bool onoff)
  { timersOn = onoff; };

  virtual int Configure(fixtures::ConfigurationType &inConfig){ 
    advancerFields = inConfig.GetValueVector<std::string>("AdvancerFields");
    //     spatialOrder   = inConfig.GetValue<int>("SpatialOrder");
    return(0); 
  };

  //  virtual int Configure(fixtures::ConfigurationType &inConfig)    { return(0); };


  virtual int InitializeAdvance()   { 
    rkH = domainPtr->TimeStep();
    return(0);
  }
  
  virtual void SetRHS(RHSType &inRHS)
  { rhsHandlerPtr = &inRHS; };

  //   pcpp::IndexIntervalType &ThreadInterval(int threadID)
  //   { return(threadIntervals[threadID]); };

  //   pcpp::IndexIntervalType &ThreadPartitionInterval(int threadID)
  //   { return(threadPartitionIntervals[threadID]); };


  void SetMessageStream(std::ostream &inStream)
  {
    messageStreamPtr = &inStream;
  };

  virtual int InitializeAdvancer(DomainType &inDomain,pcpp::ParallelGlobalType &inGlobal,std::ostream &inStream) {      

    //    FC_GLOBAL(fd1load,FD1LOAD)();

    SetRHS(inDomain.RHS());
    SetGlobal(inGlobal);
    SetMessageStream(inStream);

    domainPtr = &inDomain;
    numGrids  = inDomain.NumberOfLocalGrids();

    if(numGrids <= 0){
      inStream << "rk4advancer:Error: No grids found for domain." << std::endl;
      return(1);
    }
      
    // Grid Data
    rhsHandlers.resize(numGrids,NULL);
    gridSizes.resize(numGrids);
    numPointsGrids.resize(numGrids);
    stateMessageIds.resize(numGrids);
    gridFortranIntervals.resize(numGrids);
    gridTimes.resize(numGrids,NULL);

    inStream << "rk4advancer:Status: Found " << numGrids << " grids in domain." << std::endl;

    for(int i = 0;i < NUMTIMERS;i++){
      advancerTimers[i] = 0.0;
#ifdef USE_PAPI
      advancerCounters[i] = new long long [numCounters];
      for(int iCount = 0;iCount < numCounters;iCount++){
        advancerCounters[i][iCount] = 0;
      }
#endif
    }
    

    int numThreads = 1;

#ifdef USE_OMP
    numThreads = omp_get_max_threads();
#endif

    inStream << "rk4advancer:Status: Setting up for " << numThreads << " threads." << std::endl;

    threadStage.resize(numThreads);
    threadGrid.resize(numThreads);

    stageStates.resize(numGrids);
    stageKs.resize(numGrids);
    gridNeighbors.resize(numGrids);

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

      int iGrid = inDomain.LocalGridIndex(iLocalGrid);

      rhsHandlers[iLocalGrid] = &inDomain.RHS(iGrid);

      GridType  &domainGrid(inDomain.Grid(iGrid));
      HaloType  &gridHalo(inDomain.Grid(iGrid).Halo());
      StateType &gridState(inDomain.State(iGrid));
      StateType dummyParams;
      StateType *paramPtr = &dummyParams;

      std::vector<StateType *> &gridParamVec(inDomain.Params());
      if(gridParamVec.size() == numGrids){
        if(gridParamVec[iGrid] != NULL)
          paramPtr = gridParamVec[iGrid];
      }
      StateType &gridParams(*paramPtr);
      

      gridFortranIntervals[iLocalGrid].resize(numThreads);

      fixtures::CommunicatorType &gridComm(inDomain.Grid(iGrid).Communicator());

      int numDim = domainGrid.Dimension();
      const std::vector<double> &dX(domainGrid.DX());
      
      const std::vector<size_t> &bufferSizes(domainGrid.BufferSizes());
      size_t numPoints = domainGrid.BufferSize();


      gridSizes[iLocalGrid] = bufferSizes;
      numPointsGrids[iLocalGrid] = numPoints;

      int timeHandle = gridState.GetDataIndex("simTime");

      if(timeHandle >= 0){
        pcpp::field::dataset::DataBufferType &timeData(gridState.Field("simTime"));
        gridTimes[iLocalGrid] = timeData.Data<double>();
        inStream << "rk4advancer:Status: Using state-supplied time (" 
                 << gridTimes[iLocalGrid] << ")" << std::endl;
      } else {
        inStream << "rk4advancer:Warning: Time data not supplied in state! Using dummy time." 
                 << std::endl;
        gridTimes[iLocalGrid] = &myTime;
        *gridTimes[iLocalGrid] = 0.0;
      }
      myTime = inDomain.Time();
      *gridTimes[iLocalGrid] = myTime;
      
      rkH = domainPtr->TimeStep();
      if(rkH < 0.0){
        inStream << "rk4advancer:Error: Received negative timestep (" << rkH << ") from "
                 << "RHS object.  Aborting initialization." << std::endl;
        return(1);
      } else if (rkH == 0.0){
        inStream << "rk4advancer:Warning: Received ZERO timestep at setup. Advancer will not update time."
                 << std::endl;
      }

      numStateFields = gridState.NumStateFields();
      inStream << "rk4advancer:Status: Number of state fields = " 
               << numStateFields << std::endl;
      // Only need one intermediate state - reused for every stage
      stageStates[iLocalGrid].resize(1);
      stageKs[iLocalGrid].resize(numStages);
      stageStates[iLocalGrid][0].CopyStateData(gridState);
      for(int iStage = 0;iStage < numStages;iStage++){
        // Multiple states will be back for adjoint "revancing"
        //        stageStates[iGrid][iStage].Copy(gridState);
        stageKs[iLocalGrid][iStage].CopyStateData(gridState);
      }
      
      gridComm.CartNeighbors(gridNeighbors[iLocalGrid]);
      
      // Set up the state message
      if(communicateState){
        inStream << "rk4advancer: Setting up message for state communication."
                 << std::endl;
        stateMessageIds[iLocalGrid] = SetupStateMessage(gridState,gridHalo);
        if(stateMessageIds[iLocalGrid] < 0){
          inStream 
            << "rk4advancer:Error: Could not set up state message for Grid(" 
            << iLocalGrid+1 << "," << iGrid+1 << ")."
            << std::endl;
          return(1);
        }
      }

    }
    inStream << "rk4advancer:Status: Initialized advancer at time = " 
             << myTime << std::endl;
    return(0);    
  };
  
  int SetupStateMessage(StateType &inState,HaloType &inHalo)
  {
    int numVar = inState.NumStateVar();
    return(inHalo.CreateMessageBuffers(numVar));
  };
  
  int PackStateMessage(StateType &inState,HaloType &inHalo,
                       int messageId,int threadId)
  {
    if(threadId == 0) {
      int iVar = 0;
      for(int iField = 0;iField < numStateFields;iField++){  
        int numComponents      = inState.NumStateFieldComponents(iField);
        double *fieldData      = inState.GetStateFieldData(iField);
        inHalo.PackMessageBuffers(messageId,iVar,numComponents,fieldData);
        iVar += numComponents;
      }
    }
    return(0);
  };

  int UnpackStateMessage(StateType &inState,HaloType &inHalo,
                         int messageId,int threadId)
  {
    
    if(threadId == 0){
      int iVar = 0;
      for(int iField = 0;iField < numStateFields;iField++){  
        int numComponents      = inState.NumStateFieldComponents(iField);
        double *fieldData      = inState.GetStateFieldData(iField);
        inHalo.UnPackMessageBuffers(messageId,iVar,numComponents,fieldData);
        iVar += numComponents;
      }
    }
    return(0);
  };
  
  int SyncStates() 
  {      
    
    int myThreadId = 0;

#ifdef USE_OMP
    myThreadId = omp_get_thread_num();
#endif

    DomainType &inDomain(*domainPtr);

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

      int iGrid = inDomain.LocalGridIndex(iLocalGrid);

      GridType  &domainGrid(inDomain.Grid(iGrid));
      HaloType &gridHalo(domainGrid.Halo());
      StateType &gridState(inDomain.State(iGrid));
      fixtures::CommunicatorType &gridComm(domainGrid.Communicator());
      std::vector<int> &cartNeighbors(gridNeighbors[iLocalGrid]);
      int stateMessageId = stateMessageIds[iLocalGrid];

      if(communicateState){
        if(PackStateMessage(gridState,gridHalo,stateMessageId,myThreadId)){
          SetError("PackStateMessage",gridComm);
        }
        
        
#pragma omp barrier        
        //         if(ErrorCheck(gridComm))
        //           return(1);
#pragma omp master 
        {
          gridHalo.SendMessage(stateMessageId,cartNeighbors,gridComm);
          gridHalo.ReceiveMessage(stateMessageId,cartNeighbors,gridComm);
        }       
#pragma omp barrier
        if(UnpackStateMessage(gridState,gridHalo,stateMessageId,myThreadId)){
          SetError("UnpackStateMessage",gridComm);
        }
      }
    }
    return(0);
  };
  
  virtual int FinalizeAdvance()   
  { 
    // Update Time
#pragma omp master
    {
      myErrorCode = 0;
      myTime = myTime + rkH;
      for(int iGrid = 0;iGrid < numGrids;iGrid++)
        *gridTimes[iGrid] = *gridTimes[iGrid] + rkH;
      if(domainPtr->FinalizeAdvance(myTime))
        myErrorCode = 1;
    }
#pragma omp barrier
    return(myErrorCode);
  }
  
  virtual int GetStateUpdate(GridType &inGrid,StateType &inState,
                             StateType &stateUpdate) 
  { return(0); };
  
  virtual int InterGridExchange() { return(0); };
  
  virtual void SetError(const std::string &errorMessage,
                        fixtures::CommunicatorType &inComm){

    int myThreadId = 0;
#ifdef USE_OMP
    myThreadId = omp_get_thread_num();
#endif
    int &rkStage(threadStage[myThreadId]);

    std::ostringstream errOut;

    errOut << errorMessage << " for RKstage = " << rkStage;
    DomainAdvancerType::SetError(errOut.str(),inComm);

  };
  
  virtual int ErrorCheck(fixtures::CommunicatorType &inComm)
  { return(DomainAdvancerType::ErrorCheck(inComm)); };


  virtual int AdvanceDomain()
  {

    int myThreadId = 0;
    myErrorCode = 0;
    pcpp::ParallelGlobalType &myGlobal(*myGlobalPtr);

#ifdef USE_OMP
    myThreadId = omp_get_thread_num();
#endif


    DomainType &simulationDomain(*domainPtr);

    std::vector<GridType *>  &domainGrids(simulationDomain.Grids());
    std::vector<StateType *> &domainStates(simulationDomain.States());
    std::vector<int>         &globalGridIndices(simulationDomain.LocalGridIndices());

#pragma omp master 
    {
      myGlobal.FunctionEntry("InitializeAdvance");
    }

    if(InitializeAdvance()){
      return(error::INIT);
    }
    
#pragma omp master
    {
      myGlobal.FunctionExit("InitializeAdvance");
    }
#pragma omp barrier
    
    int &rkStage(threadStage[myThreadId]);
    int &rkGrid(threadGrid[myThreadId]);
    
    for(rkStage = 0;rkStage < numStages;rkStage++){
      
      for(rkGrid = 0;rkGrid < numGrids;rkGrid++){
        
        double tick = 0.0;
        double tock = 0.0;
        int domainGridIndex = globalGridIndices[rkGrid];
        
#ifdef USE_PAPI
        long long counterTicks[MAXCOUNTERS];
        long long counterTocks[MAXCOUNTERS];
#endif
        
        StateType &baseState(*domainStates[domainGridIndex]);
        GridType &domainGrid(*domainGrids[domainGridIndex]);
        
        fixtures::CommunicatorType &gridComm(domainGrid.Communicator());
        HaloType  &gridHalo(domainGrid.Halo());
        StateType &stageK(stageKs[rkGrid][rkStage]);
        
        // Perform a pointer jig...
        // Sets up the input state (rkInputState) and
        // the output state (rkNextState) for the current
        // RK substep.
        StateType *rkNextStatePtr = NULL;
        StateType *rkInputStatePtr = NULL;
        if(rkStage < 3) {
          rkNextStatePtr = &stageStates[rkGrid][0];
          if(rkStage == 0) {
            rkInputStatePtr = domainStates[domainGridIndex];
          } else {
            rkInputStatePtr = &stageStates[rkGrid][0];
          }
        } else { 
          rkNextStatePtr  = domainStates[domainGridIndex];
          rkInputStatePtr = &stageStates[rkGrid][0];
        } 
        StateType &rkNextState(*rkNextStatePtr);
        StateType &rkInputState(*rkInputStatePtr);
                
        std::vector<int> &cartNeighbors(gridNeighbors[rkGrid]);
        int stateMessageId = stateMessageIds[rkGrid];


#pragma omp barrier
#pragma omp master
        {
          rhsHandlerPtr = rhsHandlers[rkGrid];
          myGlobal.FunctionEntry("GetRHS");
        }
#pragma omp barrier


        //          std::cout << "INPUT STATE: " << std::endl;
        //          rkInputState.Report(std::cout);
        //          std::cout << "K: " << std::endl;
        //          stageK.Report(std::cout);
        if(GetRHS(domainGrid,gridHalo,rkInputState,stageK,myThreadId)){
          myErrorCode = 1;
        }

        //          std::cout << "RHS: " << std::endl;
        //          stageK.Report(std::cout);

#pragma omp barrier

        if(myErrorCode){
          return(error::RHS);
        }
        
#pragma omp master 
        {
          myGlobal.FunctionExit("GetRHS");
        }

        if(rkStage < 3){
#pragma omp barrier
#pragma omp master
          {
            myGlobal.FunctionEntry("ZAXPY");
          }

          AXPY(stageWeights[rkStage]*rkH,stageK,baseState,rkNextState,rkGrid,myThreadId);

#pragma omp barrier
#pragma omp master
          {
            myGlobal.FunctionExit("ZAXPY");
          }
        } else {
#pragma omp barrier
#pragma omp master
          {
            myGlobal.FunctionEntry("RKSum");
          }
          
          RKSum(rkGrid,myThreadId);
          
#pragma omp barrier
#pragma omp master
          {
            myGlobal.FunctionExit("RKSum");
          }
        }

        if(communicateState){

#pragma omp master
          {
            myGlobal.FunctionEntry("StateComm");
            myGlobal.FunctionEntry("PackMessage");
          }
          if(PackStateMessage(rkNextState,gridHalo,stateMessageId,myThreadId)){
            SetError("PackStateMessage",gridComm);
          }
#pragma omp barrier
#pragma omp master 
          {
            myGlobal.FunctionExit("PackMessage");
          }
#pragma omp barrier        
          //         if(ErrorCheck(gridComm))
          //           return(1);
#pragma omp master 
          {
            myGlobal.FunctionEntry("StateComm");
            gridHalo.SendMessage(stateMessageId,cartNeighbors,gridComm);
            gridHalo.ReceiveMessage(stateMessageId,cartNeighbors,gridComm);
            myGlobal.FunctionExit("StateComm");
          }
          
#pragma omp barrier
#pragma omp master
          {
            myGlobal.FunctionEntry("UnpackMessage");
          }

          if(UnpackStateMessage(rkNextState,gridHalo,stateMessageId,myThreadId)){
            SetError("UnpackStateMessage",gridComm);
          }
#pragma omp barrier
#pragma omp master
          {
            myGlobal.FunctionExit("UnpackMessage");
            myGlobal.FunctionExit("StateComm");
          }
          //         if(ErrorCheck(gridComm))
          //        return(1);
        }
      }
      // Intergrid exchange of substate
      if(InterGridExchange()){
        return(error::IGX); // igx does its own error handling
      //      std::cout << "rk step done" << std::endl;
      }
    }
#pragma omp master
    {
      myGlobal.FunctionEntry("FinalizeAdvance");
    }
    if(FinalizeAdvance()) 
      return(error::FINAL);
#pragma omp master
    {
      myGlobal.FunctionExit("FinalizeAdvance");
    }
    
    return(0);
    
  };
  

  
  
  int GetRHS(GridType &currentGrid,HaloType &gridHalo,StateType &inState,
             StateType &stateRHS,int myThreadId)
  {
    int &rkStage(threadStage[myThreadId]);
    
#pragma omp master
    {
      myGlobalPtr->FunctionEntry("SetupRHS");
    }

    if(rhsHandlerPtr->SetupRHS(currentGrid,inState,stateRHS,myThreadId)){
      if(messageStreamPtr)
        *messageStreamPtr << "rhsHandler::SetupRHS failed in rk stage " << rkStage
                          << std::endl;
      return(1);
    }

#pragma omp barrier
#pragma omp master
    {
      myGlobalPtr->FunctionExit("SetupRHS");
      myGlobalPtr->FunctionEntry("RHS");
    }
    int rhsErrorCode = rhsHandlerPtr->RHS(myThreadId);
    if(rhsErrorCode){
      if(messageStreamPtr)
        *messageStreamPtr << "rhsHandler::RHS failed in rk stage " << rkStage
                          << " with error code (" << rhsErrorCode << ")" << std::endl;
      return(1);
    }
    
#pragma omp barrier
#pragma omp master 
    {
      myGlobalPtr->FunctionExit("RHS");
    }
    
    return(0);
  };
  

  int AXPY(double a,StateType &xState,StateType &yState,StateType &outState,
           int localGridIndex,int myThreadId)
  {
    
#ifdef USE_PAPI
    long long counterTicks[MAXCOUNTERS];
    long long counterTocks[MAXCOUNTERS];
#endif

    DomainType &simulationDomain(*domainPtr);
    std::vector<GridType *>  &domainGrids(simulationDomain.Grids());
    std::vector<int> &globalGridIndices(simulationDomain.LocalGridIndices());
    int gridIndex = globalGridIndices[localGridIndex];
    
    GridType &opGrid(*domainGrids[gridIndex]);
    size_t numPoints = opGrid.BufferSize();

    //#if 0
    if(numPoints == 1){ // this case is needed by testing
      StateType tempState1(yState);
      StateType tempState2(xState*a);
      StateType tempState3(tempState1 + tempState2);
      outState.Copy(tempState3);
      //      outState.CopyStateData(xState*a + yState);
      //      operators::AXPY<StateType>(a,xState,outState);    
      pcpp::field::dataset::DataBufferType &timeData(outState.Field("simTime"));
      double *nextTime = timeData.Data<double>();
      myTime += a;
      *nextTime = *nextTime + a;
      return(0);
    }
    //#endif

    double tick = 0.0;
    std::vector<size_t> &myInterval(gridFortranIntervals[localGridIndex][myThreadId]);
    int numDim = opGrid.Dimension();
    std::vector<size_t> &numPointsX(opGrid.BufferSizes());

    for(int iField = 0;iField < numStateFields;iField++){

      double *xData   = xState.GetStateFieldData(iField);
      double *yData   = yState.GetStateFieldData(iField);
      double *outData = outState.GetStateFieldData(iField);

      int numComponents      = xState.NumStateFieldComponents(iField);
     
      for(int iComponent = 0;iComponent < numComponents;iComponent++){

        size_t compIndex = iComponent*numPoints;

        double *xBuffer = &xData[compIndex];
        double *yBuffer = &yData[compIndex];
        double *zBuffer = &outData[compIndex];


#pragma omp master
        {
          if(timersOn){
            tick = ix::profiler::Time();
#ifdef USE_PAPI
            PAPI_read(papiEventSet,counterTicks);
#endif
          }
        }
        
        FC_MODULE(operators,zaxpy,OPERATORS,ZAXPY)
          (&numDim,&numPoints,&numPointsX[0],&myInterval[0],
           &a,xBuffer,yBuffer,zBuffer);
        
#pragma omp master
        {
          if(timersOn){
            //        if((myThreadId == 0) && timersOn){
#ifdef USE_PAPI
            PAPI_read(papiEventSet,counterTocks);
            AccumulateCounters(counterTicks,counterTocks,advancerCounters[ZAXPYKERNEL]);
#endif           
            advancerTimers[ZAXPYKERNEL] += (ix::profiler::Time() - tick);
          }
        }
      }
    }
    return(0);
  };
  
  
  int RKSum(int localGridIndex,int myThreadId)
  {
    
#ifdef USE_PAPI
    long long counterTicks[MAXCOUNTERS];
    long long counterTocks[MAXCOUNTERS];
#endif

    DomainType &simulationDomain(*domainPtr);
    std::vector<GridType *>  &domainGrids(simulationDomain.Grids());
    std::vector<int> &globalGridIndices(simulationDomain.LocalGridIndices());
    
    int gridIndex = globalGridIndices[localGridIndex];

    GridType &opGrid(*domainGrids[gridIndex]);
    std::vector<StateType *> &domainStates(simulationDomain.States());
    StateType &baseState(*domainStates[gridIndex]);
    std::vector<StateType> &gridKs(stageKs[localGridIndex]);

    size_t numPoints = opGrid.BufferSize();
    double h = rkH;

    //#if 0
    if(numPoints == 1){
      baseState += (((gridKs[0] + gridKs[3])*(1.0/6.0) + (gridKs[1] + gridKs[2])*(1.0/3.0))*h);
      return(0);
    }
    //#endif

    double tick = 0.0;
    std::vector<size_t> &myInterval(gridFortranIntervals[localGridIndex][myThreadId]);
    std::vector<size_t> &numPointsX(opGrid.BufferSizes());
    int numDim = opGrid.Dimension();


    for(int iField = 0;iField < numStateFields;iField++){
      
      double *k1Data   = gridKs[0].GetStateFieldData(iField);
      double *k2Data   = gridKs[1].GetStateFieldData(iField);
      double *k3Data   = gridKs[2].GetStateFieldData(iField);
      double *k4Data   = gridKs[3].GetStateFieldData(iField);
      double *baseData = baseState.GetStateFieldData(iField);

      int numComponents      = baseState.NumStateFieldComponents(iField);
      
      for(int iComponent = 0;iComponent < numComponents;iComponent++){

        size_t compIndex    = iComponent*numPoints;
        double *k1Buffer    = &k1Data[compIndex];
        double *k2Buffer    = &k2Data[compIndex];
        double *k3Buffer    = &k3Data[compIndex];
        double *k4Buffer    = &k4Data[compIndex];
        double *stateBuffer = &baseData[compIndex];
        
#pragma omp master
        {
          if(timersOn){
            tick = ix::profiler::Time();
#ifdef USE_PAPI
            PAPI_read(papiEventSet,counterTicks);
#endif
          }
          
        }

        FC_MODULE(rungekutta,rk4sum,RUNGEKUTTA,RK4SUM)
          (&numDim,&numPoints,&numPointsX[0],&myInterval[0],
           &h,k1Buffer,k2Buffer,k3Buffer,k4Buffer,stateBuffer);

#pragma omp master
        {
          if(timersOn){
#ifdef USE_PAPI
            PAPI_read(papiEventSet,counterTocks);
            AccumulateCounters(counterTicks,counterTocks,advancerCounters[RKSUMKERNEL]);
#endif           
            advancerTimers[RKSUMKERNEL] += (ix::profiler::Time() - tick);
          }
        }
      }
    }
    return(0);
  };
  
  void GetTimers(double **timersPtr,const char ***namesPtr,unsigned int *numTimers)
  {
    *timersPtr  = advancerTimers;
    *namesPtr   = timerNames;
    *numTimers  = NUMTIMERS;
  };

  void GetCounters(long long ***countersPtr)
  {
    *countersPtr  = advancerCounters;
  };

  void SetGlobal(pcpp::ParallelGlobalType &inGlobal)
  {
    myGlobalPtr = &inGlobal;
  };
  
  void SyncIntervals()
  {
    
    int myThreadId = 0;
    int numThreads = 1;
    
    DomainType &simulationDomain(*domainPtr);

#ifdef USE_OMP
    myThreadId = omp_get_thread_num();
    //    numThreads = omp_get_num_threads();
#endif
    
#pragma omp barrier
    
    for(int iLocalGrid = 0;iLocalGrid < numGrids;iLocalGrid++){

      int iGrid = simulationDomain.LocalGridIndex(iLocalGrid);
      
      GridType  &inGrid(simulationDomain.Grid(iGrid));
    
      std::vector<pcpp::IndexIntervalType> &threadPartitionIntervals(inGrid.ThreadPartitionBufferIntervals());
      std::vector< std::vector<size_t> > &threadFortranIntervals(gridFortranIntervals[iGrid]);
      pcpp::IndexIntervalType &threadInterval(threadPartitionIntervals[myThreadId]);
      std::vector<size_t> &threadFortranInterval(threadFortranIntervals[myThreadId]);
      threadInterval.Flatten(threadFortranInterval);
      std::vector<size_t>::iterator tfIt = threadFortranInterval.begin();
      while(tfIt != threadFortranInterval.end()){
        *tfIt += 1;
        tfIt++;
      }
    }
    
  };
  
protected:
  
  pcpp::ParallelGlobalType *myGlobalPtr;
  bool timersOn;
  bool communicateState;

  double advancerTimers[NUMTIMERS];
  long long *advancerCounters[NUMTIMERS];

  //   bool haveBCs;
  //   bool haveHalos;

  int numStateFields;
  int myErrorCode;
  std::vector<int> threadStage;
  std::vector<int> threadGrid;

  int numGrids;
  
  int numStages;

  double rkH;

  double *timePtr;
  double *timeStepPtr;

  std::vector<RHSType *> rhsHandlers;
  RHSType *rhsHandlerPtr;

  //  int spatialOrder;
  std::vector<std::string> advancerFields;

  //  std::vector<bool>   haveBoundary;
  std::vector<double> stageWeights;
  std::vector<double> updateWeights;
  
  std::vector<std::vector<StateType> > stageStates;
  std::vector<std::vector<StateType> > stageKs;

  // Grid data
  std::vector<size_t> numPointsGrids;
  std::vector<std::vector<size_t> > gridSizes;
  std::vector<int> stateMessageIds;
  std::vector<double *> gridTimes;
  //  std::vector< std::vector<pcpp::IndexIntervalType> > gridPartitionIntervals;
  std::vector< std::vector< std::vector<size_t> > > gridFortranIntervals;

  //  StateType stateUpdate;
  
  //   std::vector<pcpp::IndexIntervalType>   threadPartitionIntervals; // local
  //   std::vector<pcpp::IndexIntervalType>   threadIntervals;          // global
  //   std::vector<std::vector<size_t> >      fortranIntervals;         // local
  
  std::vector<std::vector<int> > gridNeighbors;

  // Just in case the user does not provide data for these
  double myTime;
  double myTimeStep;
  
};