MaxwellSolver.H

Go to the documentation of this file.

#ifndef __MAXWELL_SOLVER_H__
#define __MAXWELL_SOLVER_H__

#include "Simulation.H"
#include "Stencil.H"
#include "MaxwellBC.H"



namespace Maxwell {

  int ComputeCurl(int numDim, size_t *numX, int numComponents,
                  size_t numPoints, size_t *opInterval, int *stencilID,
                  plascom2::operators::sbp::base &inOperator,
                  double *recipDeltaX, double *inputField, double *outputField);

  int ComputeMaxwellRHS_Bfield(int numDim, size_t *numX, int numComponents, size_t numPoints,
                               size_t *opInterval, int *stencilID,
                               plascom2::operators::sbp::base &inOperator,
                               double *recipDeltaX, double *Efield, double *dBdt);

  int ComputeMaxwellRHS_Dfield(int numDim, size_t *numX, int numComponents, size_t numPoints,
                               size_t *opInterval, int *stencilID,
                               plascom2::operators::sbp::base &inOperator,
                               double *recipDeltaX, double *Hfield, double *J, double *dDdt);

  int ComputeRecipEpsMu(int numDim, size_t numPoints, size_t *numX, size_t *opInterval,
                        double *epsMu, double *recipEpsMu);

  int ConvertEfieldtoDfield(int numDim, size_t numPoints, size_t *numX, size_t *opInterval,
                            double *epsMu, double *Efield, double *Dfield);

  int ConvertDfieldtoEfield(int numDim, size_t numPoints, size_t *numX, size_t *opInterval,
                            double *recipEpsMu, double *Dfield, double *Efield);

  int ConvertBfieldtoHfield(int numDim, size_t numPoints, size_t *numX, size_t *opInterval,
                            double *recipEpsMu, double *Bfield, double *Hfield);

  int ConvertHfieldtoBfield(int numDim, size_t numPoints, size_t *numX, size_t *opInterval,
                            double *epsMu, double *Hfield, double *Bfield);

  typedef simulation::grid::subregion   GridRegionType;

  template<typename GridT, typename StateT, typename OperatorT>
  class rhs : public simulation::rhs::domain_base<GridT,StateT,OperatorT> {

  public:
    typedef GridT                                                               GridType;
    typedef StateT                                                              StateType;
    typedef OperatorT                                                           OperatorType;
    typedef simulation::grid::halo                                              HaloType;
    typedef simulation::rhs::domain_base<GridType,StateType,OperatorType>       rhs_base_t;
    typedef typename rhs_base_t::BoundaryType                                   BoundaryType;
    typedef typename rhs_base_t::BCType                                         BCType;


  public:
    // - First thing to call from outside, sets up grid, halo, state and operator
    // Only called once ever
    int Initialize(GridType &inGrid, StateType &inState, const StateType &paramState, OperatorType &inOperator)
    {

      if(Init())
        return(1);

      if(SetParamState(paramState))
        return(7);

      if(SetGrid(inGrid))
        return(2);

      if(InitializeState(inState))
        return(3);

      if(HaloSetup(inGrid.Halo()))
        return(4);

      if(SetupOperator(inOperator))
        return(5);

      if(SetupTempBuffers())
        return(6);

      return(0);

    };


    // - Called any time the state/rhs are different, for RK, that's
    // every substep as each stage has separate storage
    int SetupRHS(GridType &inGrid, StateType &inState, StateType &rhsState, int myThreadID = 0)
    {

      int masterThread = (myThreadID == 0);

      if(SetGrid(inGrid))
        return(1);

      errorCode = 0;

      if(masterThread)
        errorCode = SetHalo(inGrid.Halo());

#pragma omp barrier
      if(errorCode > 0)
        return(errorCode);

      if(masterThread)
        errorCode = SetState(inState);

#pragma omp barrier
      if(errorCode > 0)
        return(errorCode);

      if(masterThread)
        errorCode = SetRHSState(rhsState);

#pragma omp barrier
      if(errorCode > 0)
        return(errorCode);

      return(0);

    };


    int HandleBoundaryConditions(int threadID)
    {

      if(!domainBCsPtr)
        return(0);
      else if (!subRegionsPtr || !domainBoundariesPtr)
        return(1);

      bool masterThread = (threadID == 0);

      std::vector<BoundaryType>   &domainBoundaries(*domainBoundariesPtr);
      std::vector<BCType>         &domainBCs(*domainBCsPtr);
      std::vector<GridRegionType> &gridRegions(*subRegionsPtr);

      int numBCs = domainBCs.size();
      int numGridRegions = gridRegions.size();
      int numDomainBoundaries = domainBoundaries.size();

      if(numBCs == 0)
        return(0);
      if(numDomainBoundaries == 0)
        return(0);

      typename std::vector<BoundaryType>::iterator dbIt = domainBoundaries.begin();

//      std::vector<double> gasParams(4,0.0);
//      gasParams[0] = cRef;
//      gasParams[1] = gamma;
//      gasParams[2] = cpRef;
//      gasParams[3] = reRef;

      while(dbIt != domainBoundaries.end()) {

        BoundaryType &domainBoundary(*dbIt++);
        GridRegionType &gridRegion(gridRegions[domainBoundary.GridRegionID()]);
        std::vector<size_t> &threadPartBufferIndices(gridRegion.threadPartitionBufferIndices[threadID]);
        size_t numBoundaryPoints = threadPartBufferIndices.size();

        if(numBoundaryPoints > 0) { // this rank/thread owns part of the boundary

          BCType &boundaryCondition(domainBoundary.BC());
          // unique int id for bc
          int bcType = boundaryCondition.BCType();

          if(bcType == simulation::domain::boundary::bc::navierstokes::HOLE)
            continue;

          int normalDir = gridRegion.normalDirection;
          std::vector<size_t> regionSizes(gridRegion.regionInterval.Sizes());
          size_t numPointsRegion = gridRegion.regionInterval.NNodes();


          // bc-associated parameters from input file
          // ///@todo get handle on BC-associated numbers/flags to avoid search
          StateType &bcParamState(boundaryCondition.Params());
          double *bcParams = bcParamState.template GetFieldBuffer<double>("Numbers");
          int *bcFlags     = bcParamState.template GetFieldBuffer<int>("Flags");

          // useful for debugging
          const std::string &bcName(boundaryCondition.BCName());
          // const std::string &regionName(gridRegion.regionName);
          // const std::string &boundaryName(domainBoundary.Name());


          // gasParams = (sndspdref,gamma,Cp,1/Re)
          // bcParams = (sigma1,sigma2,sbpBoundaryStencilWeight)
          // bcFlags = unused currently (soon to be part of BC API)

          switch(bcType) {

          case simulation::domain::boundary::bc::navierstokes::SAT_FARFIELD:
//            bcParams[2]      = boundaryWeight;
//            // Call SAT Farfield kernel
//            FC_MODULE(satutil,farfield,SATUTIL,FARFIELD)
//              ( &numDim,&bufferSizes[0],&numPointsBuffer,&normalDir,&regionSizes[0],
//                &numPointsRegion,&numBoundaryPoints,&threadPartBufferIndices[0],&gridMetric[0],
//                &gridJacobian,bcParams,&gasParams[0],myStateBuffers[0],myStateBuffers[1],
//                myStateBuffers[2],&numScalar,myStateBuffers[3],myRHSBuffers[0],myRHSBuffers[1],
//                myRHSBuffers[numDim+1],myRHSBuffers[numDim+2],targetStateBuffers[0],
//                targetStateBuffers[1],targetStateBuffers[2],targetStateBuffers[3]);
//
//            break;
            std::cout << "MaxwellRHS:HandleBoundaryConditions: Error: SAT_FARFIELD"
                      << " has not been implemented (yet)." << std::endl;
            return(1);

          case simulation::domain::boundary::bc::navierstokes::SAT_NOSLIP_ISOTHERMAL:
//            std::cout << "EulerRHS:HandleBoundaryConditions: Error: SAT_NOSLIP_ISOTHERMAL"
//                      << " has not been implemented (yet)." << std::endl;
//            return(1);
            /* Set E-field components */
            for(int iComp=0; iComp<3; iComp++) {
              double *rhsPtr = myRHSBuffers[0]+iComp*numPointsBuffer;
              double *statePtr = myStateBuffers[0]+iComp*numPointsBuffer;
//            bc::ComputeDirichletBC(threadPartBufferIndices, rhsPtr, bcParams[iComp]);
              bc::ComputeSATDirichletBC(threadPartBufferIndices, rhsPtr, statePtr, bcParams[iComp], bcParams[6]);
            }

            /* Set H-field components */
            for(int iComp=0; iComp<3; iComp++) {
              double *rhsPtr = myRHSBuffers[1]+iComp*numPointsBuffer;
              double *statePtr = myStateBuffers[3]+iComp*numPointsBuffer;
//            bc::ComputeDirichletBC(threadPartBufferIndices, rhsPtr, bcParams[iComp+3]);
              bc::ComputeSATDirichletBC(threadPartBufferIndices, rhsPtr, statePtr, bcParams[iComp+3], bcParams[7]);
            }

            break;

          case simulation::domain::boundary::bc::navierstokes::SAT_SLIP_ADIABATIC:
            std::cout << "MaxwellRHS:HandleBoundaryConditions: Error: SAT_SLIPWALL_ADIABATIC"
                      << " has not been implemented (yet)." << std::endl;
            return(1);
//            bcParams[1]      = boundaryWeight;
//            FC_MODULE(satutil,slip_adiabatic,SATUTIL,SLIP_ADIABATIC)
//              ( &numDim,&bufferSizes[0],&numPointsBuffer,&normalDir,&regionSizes[0],
//                &numPointsRegion,&numBoundaryPoints,&threadPartBufferIndices[0],&gridMetric[0],
//                &gridJacobian,bcParams,&gasParams[0],myStateBuffers[0],myStateBuffers[1],
//                myStateBuffers[2],&numScalar,myStateBuffers[3],myRHSBuffers[0],myRHSBuffers[1],
//                myRHSBuffers[numDim+1],myRHSBuffers[numDim+2],targetStateBuffers[0],
//                targetStateBuffers[1],targetStateBuffers[2],targetStateBuffers[3]);
//
//            break;

          default:
            std::cout << "ERROR! MaxwellRHS:Could not resolve BC("
                      << bcName << ") = " << bcType << std::endl;
            return(1);
            break;

          }

        }

      }

      return(0);

    };


    // No sources (yet)
    int HandleSources(int threadID)
    {
      return(0);
    };


    // main call to get RHS (SetupRHS needs to have been called first)
    int RHS(int threadID)
    {

       // Call Maxwell RHS
       if(MaxwellRHS(threadID))
        return(1);

       // Apply boundary terms
       if(HandleBoundaryConditions(threadID))
        return(2);

       // Apply additional/external source terms
       if(HandleSources(threadID))
        return(3);

       // Compute dependent variables
//       if(ComputeDV(threadID))
//      return(4);

       return(0);

    };


    // =============== Modify below at your own risk =============

    int ComputeDV(int threadID)
    {

      // Get intervals
      std::vector<size_t> &fortranBufferInterval(fortranBufferIntervals[threadID]);
      std::vector<size_t> &fortranOpInterval(fortranPartitionBufferIntervals[threadID]);

      if(ConvertHfieldtoBfield(numDim, numPointsBuffer, &bufferSizes[0], &fortranOpInterval[0],
                               myStateBuffers[5], myStateBuffers[3], myStateBuffers[1]))
      {
        std::cout << "ComputeDV - Compute B-field from H-field failed." << std::endl;
        return(1);
      }

      if(ConvertEfieldtoDfield(numDim, numPointsBuffer, &bufferSizes[0], &fortranOpInterval[0],
                               myStateBuffers[5], myStateBuffers[0], myStateBuffers[2]))
      {
        std::cout << "ComputeDV - Compute D-field from E-field failed." << std::endl;
        return(2);
      }

      return(0);

    };


    int MaxwellRHS(int threadID)
    {

      // Get intervals
      std::vector<size_t> &fortranBufferInterval(fortranBufferIntervals[threadID]);
      std::vector<size_t> &fortranOpInterval(fortranPartitionBufferIntervals[threadID]);

      //       std::cout << "Buffer Sizes: (";
      //       pcpp::io::DumpContents(std::cout,bufferSizes,",");
      //       std::cout << ")" << std::endl;
      //       std::cout << "Fortran Buffer Interval: (";
      //       pcpp::io::DumpContents(std::cout,fortranBufferInterval,",");
      //       std::cout << ")" << std::endl;
      //       std::cout << "Fortran OpInterval: (";
      //       pcpp::io::DumpContents(std::cout,fortranOpInterval,",");
      //       std::cout << ")" << std::endl;
      //       std::cout << "BField(before):" << std::endl;
      //       for(size_t iPoint = 0;iPoint < numDim*numPointsBuffer;iPoint++)
      //         std::cout << bFieldPtr[iPoint] << " ";
      //       std::cout << std::endl;
      //       std::cout << "dBdt(BEFORE):" << std::endl;
      //       for(size_t iPoint = 0;iPoint < numDim*numPointsBuffer;iPoint++)
      //         std::cout << myRHSBuffers[1][iPoint] << " ";
      //       std::cout << std::endl;
      //       size_t iPoint = 0;
      //       for(int iDim = 0;iDim < numDim;iDim++){
      //         std::cout << "StencilID[" << iDim << "]:" << std::endl;
      //         for(size_t iZ = 0;iZ < bufferSizes[2];iZ++){
      //           for(size_t iY = 0;iY < bufferSizes[1];iY++){
      //             for(size_t iX = 0;iX < bufferSizes[0];iX++){
      //               std::cout << stencilID[iPoint++] << " ";
      //             }
      //           }
      //         }
      //       }

      // Compute RHS for B-field
      if(ComputeMaxwellRHS_Bfield(numDim, &bufferSizes[0], numComponents, numPointsBuffer,
                                  &fortranOpInterval[0], stencilID, myOperator, &recipDX[0],
                                  myStateBuffers[0], myRHSBuffers[1]))
      {
        return(1);
      }

      //       std::cout << "dBdt(AFTER):" << std::endl;
      //       for(size_t iPoint = 0;iPoint < numDim*numPointsBuffer;iPoint++)
      //         std::cout << myRHSBuffers[1][iPoint] << " ";
      //       std::cout << std::endl;

      // Convert dB/dt to dH/dt
      if(ConvertBfieldtoHfield(numDim, numPointsBuffer, &bufferSizes[0], &fortranOpInterval[0],
                               &recipEpsMu[0], myRHSBuffers[1], myRHSBuffers[1]))
      {
        return(2);
      }

      //       std::cout << "dHdt:" << std::endl;
      //       for(size_t iPoint = 0;iPoint < numDim*numPointsBuffer;iPoint++)
      //         std::cout << myRHSBuffers[1][iPoint] << " ";
      //       std::cout << std::endl;
      
      // Compute RHS for D-field
      if(ComputeMaxwellRHS_Dfield(numDim, &bufferSizes[0], numComponents, numPointsBuffer,
                                  &fortranOpInterval[0], stencilID, myOperator, &recipDX[0],
                                  myStateBuffers[3], myStateBuffers[4], myRHSBuffers[0]))
      {
        return(3);
      }

      // Convert dD/dt to dE/dt
      if(ConvertDfieldtoEfield(numDim, numPointsBuffer, &bufferSizes[0], &fortranOpInterval[0],
                               &recipEpsMu[0], myRHSBuffers[0], myRHSBuffers[0]))
      {
        return(4);
      }

      return(0);

    };


    std::vector<double *> &StateBuffers()
    {
      return(myStateBuffers);
    };


    std::vector<double *> &RHSBuffers()
    {
      return(myRHSBuffers);
    };


    int SetParamState(const StateType &paramState)
    {

      int myThreadID = 0;
      int numThreads = 1;
      int masterThread = 1;

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

      if(masterThread) {
        epsilon0Handle = paramState.GetDataIndex("epsilon0");
        if(epsilon0Handle < 0)
          return(1);
        const pcpp::field::dataset::DataBufferType &epsilon0Data(paramState.Field(epsilon0Handle));
        epsilon0Ptr = epsilon0Data.Data<double>();
        if(!epsilon0Ptr)
          return(1);
        epsilon0 = *epsilon0Ptr;

        mu0Handle = paramState.GetDataIndex("mu0");
        if(mu0Handle < 0)
          return(2);
        const pcpp::field::dataset::DataBufferType &mu0Data(paramState.Field(mu0Handle));
        mu0Ptr = mu0Data.Data<double>();
        if(!mu0Ptr)
          return(2);
        mu0 = *mu0Ptr;

        cflHandle = paramState.GetDataIndex("inputCFL");
        if(cflHandle < 0)
          return(3);
        const pcpp::field::dataset::DataBufferType &cflData(paramState.Field(cflHandle));
        cflPtr = cflData.Data<double>();
        if(!cflPtr)
          return(3);
        cfl = *cflPtr;

        deltaTHandle = paramState.GetDataIndex("inputDT");
        if(deltaTHandle < 0)
          return(4);
        const pcpp::field::dataset::DataBufferType &deltaTData(paramState.Field(deltaTHandle));
        deltaTPtr = deltaTData.Data<double>();
        if(!deltaTPtr)
          return(4);
        deltaT = *deltaTPtr;

        cRefHandle = paramState.GetDataIndex("refC");
        if(cRefHandle < 0)
          return(5);
        const pcpp::field::dataset::DataBufferType &cRefData(paramState.Field(cRefHandle));
        cRefPtr = cRefData.Data<double>();
        if(!cRefPtr)
          return(5);
        cRef = *cRefPtr;
      }
#pragma omp barrier

      return(0);

    };


    int InitializeState(StateType &inState)
    {

      if(!gridInitialized)
        return(1);

      int myThreadID = 0;
      int numThreads = 1;
      int masterThread = 1;

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

      if(masterThread) {
        timeHandle = inState.GetDataIndex("simTime");
        if(timeHandle < 0)
          return(1);
        pcpp::field::dataset::DataBufferType &timeData(inState.Field(timeHandle));
        timePtr = timeData.Data<double>();
        if(!timePtr)
          return(1);
        time = *timePtr;

        eFieldHandle = inState.GetFieldHandle("Efield");
        if(eFieldHandle < 0)
          return(1);
        pcpp::field::dataset::DataBufferType &eFieldData(inState.GetFieldDataByHandle(eFieldHandle));
        eFieldPtr    = eFieldData.Data<double>();
        if(!eFieldPtr)
          return(1);
    
        bFieldHandle = inState.GetFieldHandle("Bfield");
        if(bFieldHandle < 0)
          return(1);
        pcpp::field::dataset::DataBufferType &bFieldData(inState.GetFieldDataByHandle(bFieldHandle));
        bFieldPtr   = bFieldData.Data<double>();
        if(!bFieldPtr)
          return(1);

        dFieldHandle = inState.GetDataIndex("Dfield");
        if(dFieldHandle < 0)
          return(1);
        pcpp::field::dataset::DataBufferType &dFieldData(inState.Field(dFieldHandle));
        dFieldPtr   = dFieldData.Data<double>();
        if(!dFieldPtr)
          return(1);

        hFieldHandle = inState.GetDataIndex("Hfield");
        if(hFieldHandle < 0)
          return(1);
        pcpp::field::dataset::DataBufferType &hFieldData(inState.Field(hFieldHandle));
        hFieldPtr   = hFieldData.Data<double>();
        if(!hFieldPtr)
          return(1);

        jHandle = inState.GetDataIndex("J");
        if(jHandle < 0)
          return(1);
        pcpp::field::dataset::DataBufferType &jData(inState.Field(jHandle));
        jPtr   = jData.Data<double>();
        if(!jPtr)
          return(1);

        epsilonMuHandle = inState.GetDataIndex("epsilonMu");
        if(epsilonMuHandle < 0)
          return(1);
        pcpp::field::dataset::DataBufferType &epsilonMuData(inState.Field(epsilonMuHandle));
        epsilonMuPtr    = epsilonMuData.Data<double>();
        if(!epsilonMuPtr)
          return(1);
#if 0
      epsilonHandle = inState.GetDataIndex("epsilon");
      if(epsilonHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &epsilonData(inState.Field(epsilonHandle));
      epsilonPtr    = epsilonData.Data<double>();
      if(!epsilonPtr)
        return(1);

      muHandle = inState.GetDataIndex("mu");
      if(muHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &muData(inState.Field(muHandle));
      muPtr    = muData.Data<double>();
      if(!muPtr)
        return(1);
#endif
        myStateBuffers.resize(6,NULL);
        myStateBuffers[0] = eFieldPtr;
        myStateBuffers[1] = bFieldPtr;
        myStateBuffers[2] = dFieldPtr;
        myStateBuffers[3] = hFieldPtr;
        myStateBuffers[4] = jPtr;
        myStateBuffers[5] = epsilonMuPtr;
//      myStateBuffers[3] = epsilonPtr;
//      myStateBuffers[4] = muPtr;
        myRHSBuffers.resize(2,NULL);

        stateInitialized = true;
      }
#pragma omp barrier

      return(0);

    };


    int SetupStencilConnectivities()
    {

      if(!operatorInitialized)
        return(1);

      int myThreadID = 0;
      int numThreads = 1;
      int masterThread = 1;

#ifdef USE_OMP
      myThreadID = omp_get_thread_num();
      numThreads = omp_get_num_threads();
      masterThread = (myThreadID == 0);
#endif

      if(masterThread) {
        stencilID = new int [numDim*numPointsBuffer];
      }

#pragma omp barrier

      // Create stencil connectivity
      if(numThreads == 1) {
        if(plascom2::operators::sbp::CreateStencilConnectivity(numDim,&bufferSizes[0],&opInterval[0],boundaryDepth,
                                                               &periodicDirs[0],stencilID,true))
              {
          std::cout << "CreateStencilConnectivity failed." << std::endl;
          return(1);
              }
      } else {
        if(plascom2::operators::sbp::CreateStencilConnectivity(numDim,&bufferSizes[0],&fortranPartitionBufferIntervals[myThreadID][0],
                                                               boundaryDepth,&periodicDirs[0],stencilID,true))
              {
          std::cout << "CreateStencilConnectivity failed." << std::endl;
          return(1);
              }
      }


      // plascom2::operators::sbp::CreateStencilConnectivity(numDim,&gridSize[0],opInterval,boundaryDepth,stencilID,true);
      
      // bug! unlike stencilID, these need to be thread-aware
      // @todo Thread-aware dual stencil connectivities in RHS
      // Invert the connectivity
      // plascom2::operators::InvertStencilConnectivity(numDim,opInterval,stencilID,dualStencilConn,true);
      // dualStencilConn  = new size_t [numDim*numPointsBuffer];        // [DIM1[stencil1Points,stencil2Points....stencil(numStencils)Points]...]
      // numPointsStencil = new size_t [numDim*numStencils]; // [DIM1[numPoints1,numPoints2....numPoints(numStencils)]...]
     
//      if(plascom2::operators::sbp::InvertStencilConnectivity(numDim,&bufferSizes[0],&opInterval[0],numStencils,
//                                                             stencilID,dualStencilConn,numPointsStencil,true))
//      {
//        std::cout << "InvertStencilConnectivity failed." << std::endl;
//        return(1);
//      }

      return(0);

    };


    int SetupOperator(OperatorType &inOperator)
    {

      int myThreadID = 0;
      int numThreads = 1;
      int masterThread = 1;

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

      if(masterThread && !operatorInitialized) {
        myOperator.Copy(inOperator);

        //      plascom2::operators::sbp::Initialize(myOperator,interiorSpatialOrder);
        numStencils    = myOperator.numStencils;
        stencilStarts  = myOperator.stencilStarts;
        stencilOffsets = myOperator.stencilOffsets;
        stencilSizes   = myOperator.stencilSizes;
        stencilWeights = myOperator.stencilWeights;
        boundaryDepth  = myOperator.boundaryDepth; // (numStencils - 1)/2;
        boundaryWeight = myOperator.boundaryWeight;
        numComponents  = 1;
        //         std::cout << "NumStencils: " << numStencils << std::endl
        //                   << "StencilStarts/sizes: (";
        //         for(int iStencil = 0;iStencil < numStencils;iStencil++){
        //           std::cout << stencilStarts[iStencil] << "/" << stencilSizes[iStencil];
        //           if(iStencil != (numStencils - 1)) std::cout << ",";
        //         }
        //         std::cout << ") " << std::endl;
        // Forticate the stencil starts
        for(int iStencil=0; iStencil<numStencils; iStencil++)
          stencilStarts[iStencil]++;

        operatorInitialized = true;
      }
#pragma omp barrier

      if(SetupStencilConnectivities())
        return(1);

      return(0);
      
    };

    int SetState(StateType &inState)
    {

      if(!stateInitialized)
        return(1);

      eFieldPtr         = inState.GetRealFieldBufferByHandle(eFieldHandle);
//      bFieldPtr         = inState.GetRealFieldBuffer(bFieldHandle);
//      dFieldPtr         = inState.GetRealFieldBuffer(dFieldHandle);
      hFieldPtr         = inState.GetRealFieldBufferByHandle(hFieldHandle);
//      jPtr              = inState.GetRealFieldBuffer(jHandle);
//      epsilonMuPtr    = inState.GetRealFieldBuffer(epsilonMuHandle);
//      epsilonPtr      = inState.GetRealFieldBuffer(epsilonHandle);
//      muPtr           = inState.GetRealFieldBuffer(muHandle);

      myStateBuffers[0] = eFieldPtr;
//      myStateBuffers[1] = bFieldPtr;
//      myStateBuffers[2] = dFieldPtr;
      myStateBuffers[3] = hFieldPtr;
//      myStateBuffers[4] = jPtr;
//      myStateBuffers[5] = epsilonMuPtr;
//      myStateBuffers[3] = epsilonPtr;
//      myStateBuffers[4] = muPtr;

      // Compute reciprocals of epsilon and mu
      if(dynamicEpsMu || recipEpsMu.empty()) {
        recipEpsMu.resize(2*numPointsBuffer);
        if(ComputeRecipEpsMu(numDim, numPointsBuffer, &bufferSizes[0], &opInterval[0],
                             myStateBuffers[5], &recipEpsMu[0]))
        {
            return(1);
        }
      }

#pragma omp barrier

      return(0);

    };


    int InitThreadIntervals()
    {

      int myThreadId = 0;
      int numThreads = 1;
      int masterThread = 1;

#ifdef USE_OMP
      myThreadId = omp_get_thread_num();
      numThreads = omp_get_max_threads();
      masterThread = (myThreadId == 0);
#endif
      
      if(masterThread) {
        fortranBufferIntervals.resize(numThreads);
        fortranPartitionBufferIntervals.resize(numThreads);
      }
      
#pragma omp barrier
      
      std::vector<pcpp::IndexIntervalType> &threadPartitionBufferIntervals(gridPtr->ThreadPartitionBufferIntervals());
      if(threadPartitionBufferIntervals.size() != numThreads)
        return(1);
      threadPartitionBufferIntervals[myThreadId].Flatten(fortranPartitionBufferIntervals[myThreadId]);

      std::vector<pcpp::IndexIntervalType> &threadBufferIntervals(gridPtr->ThreadBufferIntervals());
      if(threadBufferIntervals.size() != numThreads)
        return(1);
      threadBufferIntervals[myThreadId].Flatten(fortranBufferIntervals[myThreadId]);
      
      // Forticate the thread intervals
      std::vector<size_t>::iterator opIt = fortranPartitionBufferIntervals[myThreadId].begin();
      while(opIt != fortranPartitionBufferIntervals[myThreadId].end()) {
        *opIt += 1;
        opIt++;
      }

      std::vector<size_t>::iterator bufIt = fortranBufferIntervals[myThreadId].begin();
      while(bufIt != fortranBufferIntervals[myThreadId].end()) {
        *bufIt += 1;
        bufIt++;
      }
      
      return(0);

    };


    int SetGrid(GridType &inGrid)
    {

      int numThreads = 1;
#ifdef USE_OMP
        numThreads = omp_get_num_threads();
#endif

      if(gridPtr != &inGrid) {
        int myThreadID = 0;
        int masterThread = 1;

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

        if(masterThread) {
          gridPtr                 = &inGrid;
          gridSizes               = gridPtr->GridSizes();
          bufferSizes             = gridPtr->BufferSizes();
          partitionInterval       = gridPtr->PartitionInterval();
          partitionBufferInterval = gridPtr->PartitionBufferInterval();
          numPointsPartition      = partitionInterval.NNodes();
          numDim                  = gridSizes.size();
          numPointsBuffer         = gridPtr->BufferSize();
          gridCommPtr             = &gridPtr->Communicator();
          subRegionsPtr           = &gridPtr->SubRegions();
          gridCommPtr->CartNeighbors(gridNeighbors);
          gridInitialized         = true;
          if(partitionBufferInterval.empty()) {
            return(1);
          }
          partitionBufferInterval.Flatten(opInterval);

          // Get dX
          dX = gridPtr->DX();

          // Allocate memory for reciprocal of dX
          recipDX.resize(numDim*numPointsBuffer);

          // Compute reciprocal of dX
          if(ComputeRecipDXUniformGrid())
          {
            return(2);
          }

          // Forticate the opInterval
          std::vector<size_t>::iterator opIt = opInterval.begin();
          while(opIt != opInterval.end()) {
            *opIt += 1;
            opIt++;
          }
        }

        // Reset the thread intervals for this grid
        InitThreadIntervals();
      }

      return(0);

    };


    int SetRHSState(StateType &rhsState)
    {

      if(dEHandle < 0) {
        dEHandle  = rhsState.GetFieldHandle("Efield");
        dHHandle  = rhsState.GetFieldHandle("Hfield");
        //        dJHandle  = rhsState.GetDataIndex("J");
      }

      eFieldRHSPtr         = rhsState.GetRealFieldBufferByHandle(dEHandle);
      hFieldRHSPtr         = rhsState.GetRealFieldBufferByHandle(dHHandle);

      myRHSBuffers[0] = eFieldRHSPtr;
      myRHSBuffers[1] = hFieldRHSPtr;
#if 0
      myRHSBuffers[0] = eFieldRHSPtr;
      myRHSBuffers[1] = &eFieldRHSPtr[numPoints];
      myRHSBuffers[2] = &eFieldRHSPtr[2*numPoints];
      myRHSBuffers[3] = bFieldRHSPtr;
      myRHSBuffers[4] = &bFieldRHSPtr[numPoints];
      myRHSBuffers[5] = &bFieldRHSPtr[2*numPoints];
#endif
      return(0);    

    };


    int HaloSetup(HaloType &inHalo)
    {

      int myThreadID = 0;
      int numThreads = 1;
      int masterThread = 1;

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

      if(masterThread) {
        SetHalo(inHalo);
        periodicDirs = inHalo.PeriodicDirs();
        if(periodicDirs.empty()) {
          std::cout << "WARNING: Defaulting to periodic domain in RHS object." << std::endl;
          periodicDirs.resize(numDim,1);
        }
      }
#pragma omp barrier

      return(0);

    };


    int SetHalo(HaloType &inHalo)
    {

      haloPtr = &inHalo;

      return(0);

    };


    // This function sets the domainBoundaries data structure *and* goes
    // through the stencilConnectivity to assign boundary stencils.  It is
    // redundant for the actual grid/domain boundaries, but boundary stencils
    // for internal boundaries are set here.
    void SetDomainBoundaries(std::vector<BoundaryType> &domainBoundaries)
    {

      domainBoundariesPtr = &domainBoundaries;

      assert(subRegionsPtr != NULL);

      std::vector<GridRegionType> &gridRegions(*subRegionsPtr);
      int numDomainBoundaries = domainBoundaries.size();
      int numGridRegions      = gridRegions.size();
      int holeStencil         = myOperator.numStencils;

      typename std::vector<BoundaryType>::iterator dbIt = domainBoundaries.begin();

      while(dbIt != domainBoundaries.end()) {

        BoundaryType &domainBoundary(*dbIt++);
        GridRegionType &gridRegion(gridRegions[domainBoundary.GridRegionID()]);

        int boundaryNormalDirection = gridRegion.normalDirection;
        pcpp::IndexIntervalType &boundaryPartitionInterval(gridRegion.regionPartitionInterval);

        int returnCode = 0;
         if(boundaryNormalDirection != 0) {
           returnCode = plascom2::operators::sbp::BoundaryStencilConnectivity(bufferSizes,
                                                                              partitionInterval,
                                                                              partitionBufferInterval,
                                                                              boundaryPartitionInterval,
                                                                              boundaryDepth,
                                                                              boundaryNormalDirection,
                                                                              stencilID);
         } else {
           returnCode = plascom2::operators::sbp::HoleStencilConnectivity(bufferSizes,
                                                                          partitionInterval,
                                                                          partitionBufferInterval,
                                                                          boundaryPartitionInterval,
                                                                          holeStencil,
                                                                          stencilID);
         }

        assert(returnCode == 0);

      }

    };


    void SetDomainBCs(std::vector<BCType> &domainBCs)
    { domainBCsPtr = &domainBCs; };


    int InitStep(GridType &inGrid,StateType &inState)
    {
      return(0);
    };


    int SetupTempBuffers()
    {
      return(0);
    };


    int ComputeRecipDXUniformGrid()
    {

      size_t iDimxnumPointsBuffer;
      double oneOverDX;

      if (numDim != 3) return(1);

      for (int iDim=0; iDim<numDim; iDim++) {
        iDimxnumPointsBuffer = iDim*numPointsBuffer;
        oneOverDX = 1.0/dX[iDim];
        for (size_t iPoint=0; iPoint<numPointsBuffer; iPoint++) {
          recipDX[iPoint+iDimxnumPointsBuffer] = oneOverDX;
        }
      }

      return(0);

    };
    double TimeStep(double *outCFL = NULL)
    {
      return(deltaT);
    };

  private:
    int Init()
    {

      int myThreadID = 0;
      int numThreads = 1;
      int masterThread = 1;

#ifdef USE_OMP
      myThreadID = omp_get_thread_num();
      numThreads = omp_get_num_threads();
      masterThread = (myThreadID == 0);
#endif

      if(masterThread) {
        // Dynamic material properties
        dynamicEpsMu = false;

        // Conserved state
        eFieldHandle  = -1;
        bFieldHandle  = -1;
        dFieldHandle  = -1;
        hFieldHandle  = -1;
        jHandle       = -1;
        dEHandle      = -1;
        dHHandle      = -1;

        // Parameters
        timeHandle          = -1;
        epsilonMuHandle          = -1;
        gridSpacingHandle = -1;
        epsilon0Handle    = -1;
        mu0Handle         = -1;
        cflHandle         = -1;
        deltaTHandle      = -1;
        cRefHandle        = -1;

        // Misc
        numPointsPartition = 0;
        numPointsBuffer    = 0;
        numDim             = 0;
        numStateFields     = 3;
        numStateVar        = 9;

        // Pointers
        gridPtr             = NULL;
        timePtr             = NULL;
        eFieldPtr           = NULL;
        bFieldPtr           = NULL;
        dFieldPtr           = NULL;
        hFieldPtr           = NULL;
        jPtr                = NULL;
        epsilonMuPtr        = NULL;
        eFieldRHSPtr        = NULL;
        hFieldRHSPtr        = NULL;
        rhsHaloPtr          = NULL;
        epsilon0Ptr         = NULL;
        mu0Ptr              = NULL;
        cflPtr              = NULL;
        deltaTPtr           = NULL;
        cRefPtr             = NULL;
        domainBoundariesPtr = NULL;
        subRegionsPtr       = NULL;
        domainBCsPtr        = NULL;

        stateInitialized    = false;
        gridInitialized     = false;
        operatorInitialized = false;
        errorCode = 0;
      }

      return(0);

    };


  protected:
    // Data handles
    int eFieldHandle;
    int bFieldHandle;
    int dFieldHandle;
    int hFieldHandle;
    int jHandle;
    int dEHandle;
    int dHHandle;
    int timeHandle;
    int epsilonMuHandle;
    int gridSpacingHandle;
    int epsilon0Handle;
    int mu0Handle;
    int cflHandle;
    int deltaTHandle;
    int cRefHandle;

    // Parameters
    double time;
    double epsilon0;
    double mu0;
    double cfl;
    double deltaT;
    double cRef;
    int numStateFields;
    int numStateVar;

    // Grid data
    GridType *gridPtr;
    fixtures::CommunicatorType *gridCommPtr;
    std::vector<int> gridNeighbors;
    int numDim;
    std::vector<size_t> gridSizes;
    std::vector<size_t> bufferSizes;
    pcpp::IndexIntervalType partitionInterval;
    pcpp::IndexIntervalType partitionBufferInterval;
    size_t numPointsPartition;
    size_t numPointsBuffer;
    std::vector<double> gridMetric;
    double gridJacobian;
    std::vector<double> dX;
    std::vector<int> periodicDirs;
    std::vector<size_t> opInterval;
    HaloType *haloPtr;
    std::vector<std::vector<size_t> > fortranBufferIntervals;
    std::vector<std::vector<size_t> > fortranPartitionBufferIntervals;

    // Parameter information
    const double *epsilon0Ptr;
    const double *mu0Ptr;
    const double *cflPtr;
    const double *deltaTPtr;
    const double *cRefPtr;

    // Input state information
    double *eFieldPtr;
    double *bFieldPtr;
    double *dFieldPtr;
    double *hFieldPtr;
    double *jPtr;
    double *epsilonMuPtr;
    double *timePtr;

    // RHS state information
    //  double *rhsTimePtr;
    double *eFieldRHSPtr;
    double *hFieldRHSPtr;

    // My storage
    std::vector<double *> myStateBuffers;
    std::vector<double *> myRHSBuffers;

    // Operator
    OperatorType myOperator;
    int numStencils;
    int boundaryDepth;
    int numStencilValues;
    int numComponents;
    int *stencilOffsets;
    int *stencilID;
    int *stencilStarts;
    int *stencilSizes;
    double *stencilWeights;
    size_t *dualStencilConn;
    size_t *numPointsStencil;
//    std::vector<size_t> opInterval; // forticated, flat partition buffer interval

    // Bookkeeping
    bool stateInitialized;
    bool gridInitialized;
    bool operatorInitialized;
    std::vector<bool> haveHalo;
    int errorCode;

    // Halo
    HaloType *rhsHaloPtr;

    // Boundary condition support
    std::vector<BoundaryType> *domainBoundariesPtr;
    std::vector<GridRegionType> *subRegionsPtr;
    std::vector<BCType> *domainBCsPtr;
    double boundaryWeight;

    // Temporary buffers
    std::vector<double> recipEpsMu;
    std::vector<double> recipDX;
    bool dynamicEpsMu;

  };

}
#endif