EulerRHS.H

Go to the documentation of this file.

#ifndef __EULER_RHS_H__
#define __EULER_RHS_H__

#include <limits>
#ifdef USE_OMP
#include <omp.h>
#endif

#include "Simulation.H"
#include "EulerUtil.H"
#include "EulerKernels.H"
#include "ViscidUtil.H"
#include "ViscidKernels.H"
#include "OperatorKernels.H"
#include "MetricKernels.H"
#include "GridKernels.H"
#include "Stencil.H"
#include "SATKernels.H"
#include "WENO.H"
#include "EOS.H"

//static const int ALTERNATERHS = 1;
static const int EXCHANGEFLUX = 2;
static const int EXCHANGEDV   = 4;
static const int ARTDISS      = 8;
static const int USEWENO      = 16;

static const int holeMask     = 1 << mask::HOLE;

using namespace simulation::domain::boundary::bc;
#define CONSTANT_DT  0
#define CONSTANT_CFL 1
namespace error {
  static const int DV       = 1;
  static const int TV       = 2;
  static const int VELGRAD  = 4;
  static const int VELDIV   = 8;
  static const int TGRAD    = 16;
  static const int CONCGRAD = 32;
  static const int TAU      = 64;
  static const int HEATFLUX = 128;
  static const int AO       = 256;
  static const int XANGE    = 512;
  static const int IFLUX    = 1024;
  static const int VFLUX    = 2048;
  static const int VEFLUX   = 4096;
  static const int BC       = 8192;
}

namespace euler {

  enum GASMODELS {IDEAL=0,NUM_GAS_MODELS};
  enum TRANSPORTMODELS {POWER=0,NUM_TRANSPORT_MODELS};

  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;
    typedef typename rhs_base_t::GlobalType   GlobalType;

    using rhs_base_t::globalPtr;

  public:

    rhs() : rhs_base_t(), exchangeFlux(true), 
            exchangeDV(false), artificialDissipation(false),useWENO(false)
    {
      globalPtr = NULL;
    };


    // ===== Main interface functions ====

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

      tauCounter = 0;

      if(Init()){
        std::cout << "EulerRHS:Fatal error: Init failed." << std::endl;
        return(1);
      }

      if(SetGrid(inGrid)){
        std::cout << "EulerRHS:Fatal error: Set Grid failed." << std::endl;
        return(2);
      }

      if(SetParamState(inParam)){
        std::cout << "EulerRHS:Fatal error: Set params failed." << std::endl;
        return(6);
      }
      
      if(globalPtr){
#pragma omp master
        {
          globalPtr->FunctionEntry("InitializeState");
        }
      }
      if(InitializeState(inState)){
        std::cout << "EulerRHS:Fatal error: Initialize state failed." << std::endl;
        return(3);
      }


      if(SetTargetState(inState)){
        if(globalPtr){
          globalPtr->ErrOut("Target data setup failed.\n");
        } else {
          std::cout << "EulerRHS::InitState:Error: Target data setup failed."
                    << std::endl;
        }
        return(1);
      }

      if(globalPtr){
#pragma omp master
        {
          globalPtr->FunctionExit("InitializeState");
          //          globalPtr->FunctionEntry("InitMinSpacing");
        }
      }

      //      InitMinSpacing();

      if(globalPtr){
#pragma omp master
        {
          //          globalPtr->FunctionExit("InitMinSpacing");
          globalPtr->FunctionEntry("HaloSetup");
        }
      }

      if(HaloSetup(inGrid.Halo())){
        std::cout << "EulerRHS:Fatal error: Halo setup failed." << std::endl;
        return(4);
      }

      if(globalPtr){
#pragma omp master
        {
          globalPtr->FunctionExit("HaloSetup");
          globalPtr->FunctionEntry("SetupOperators");
        }
      }

      if(SetupOperators(inOperator)){
        std::cout << "EulerRHS:Fatal error: Operator setup failed." << std::endl;
        return(5);
      }
      if(globalPtr){
#pragma omp master
        {
          globalPtr->FunctionExit("SetupOperators");
        }
      }


      //      std::cout << "Configured with grid type = " << gridType << std::endl;
      return(0);

    };

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

      tauCounter = 0;

      if(Init()){
        std::cout << "EulerRHS:Fatal error: Init failed." << std::endl;
        return(1);
      }

      if(SetGrid(inGrid)){
        std::cout << "EulerRHS:Fatal error: Set Grid failed." << std::endl;
        return(2);
      }

      if(SetParamState(inParam)){
        std::cout << "EulerRHS:Fatal error: Set params failed." << std::endl;
        return(6);
      }
      

      if(globalPtr){
#pragma omp master
        {
          globalPtr->FunctionEntry("InitializeState");
        }
      }
      if(InitializeState(inState)){
        std::cout << "EulerRHS:Fatal error: Initialize state failed." << std::endl;
        return(3);
      }

      if(SetTargetState(targetState)){
        if(globalPtr){
          globalPtr->ErrOut("Target data setup failed.\n");
        } else {
          std::cout << "EulerRHS::InitState:Error: Target data setup failed."
                    << std::endl;
        }
        return(1);
      }

      if(globalPtr){
#pragma omp master
        {
          globalPtr->FunctionExit("InitializeState");
          //          globalPtr->FunctionEntry("InitMinSpacing");
        }
      }

      //      InitMinSpacing();

      if(globalPtr){
#pragma omp master
        {
          //          globalPtr->FunctionExit("InitMinSpacing");
          globalPtr->FunctionEntry("HaloSetup");
        }
      }

      if(HaloSetup(inGrid.Halo())){
        std::cout << "EulerRHS:Fatal error: Halo setup failed." << std::endl;
        return(4);
      }

      if(globalPtr){
#pragma omp master
        {
          globalPtr->FunctionExit("HaloSetup");
          globalPtr->FunctionEntry("SetupOperators");
        }
      }

      if(SetupOperators(inOperator)){
        std::cout << "EulerRHS:Fatal error: Operator setup failed." << std::endl;
        return(5);
      }
      if(globalPtr){
#pragma omp master
        {
          globalPtr->FunctionExit("SetupOperators");
        }
      }

      //      std::cout << "Configured with grid type = " << gridType << std::endl;
      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)
    {
      
      
      errorCode = 0;

#pragma omp master
      {
        errorCode = SetHalo(inGrid.Halo());
      }

#pragma omp barrier    
      
      if(errorCode > 0)
        return(errorCode);
      
#pragma omp master
      {
        errorCode = SetState(inState);
      }      
#pragma omp barrier    
      
      if(errorCode > 0)
        return(errorCode);
      
#pragma omp master
      {
        errorCode = SetRHSState(rhsState);
      }      

#pragma omp barrier    

      if(errorCode > 0)
        return(errorCode);
      
      return(0);
      
    };
    
    // main call to get RHS (SetupRHS needs to have been called first)
    int PrepareBuffers(int threadId)                                              
    {

      {
        std::vector<size_t> 
          &fortranBufferInterval(fortranBufferIntervals[threadId]);
        
        // Zero the RHS buffers
        if(globalPtr){
#pragma omp barrier
#pragma omp master
          {
            globalPtr->FunctionEntry("ZeroRHS");
          }
        }

        double zero = 0.0;
        for(int iEquation = 0;iEquation < numEquations;iEquation++){
          FC_MODULE(operators,assignmentxa,OPERATORS,ASSIGNMENTXA)
            (&numDim,&numPointsBuffer,&bufferSizes[0],
             &fortranBufferInterval[0],&zero,myRHSBuffers[iEquation]);

          if(iEquation == 0){
            FC_MODULE(operators,assignmentxa,OPERATORS,ASSIGNMENTXA)
              (&numDim,&numPointsBuffer,&bufferSizes[0],
               &fortranBufferInterval[0],&zero,viscidFluxBuffers[iEquation]);
          }

        }

        if(globalPtr){
#pragma omp barrier
#pragma omp master
          {
            globalPtr->FunctionExit("ZeroRHS");
          }
        }
      }
      return(0);
    };
    
    int RHS(int threadId)
    {

      PrepareBuffers(threadId);

      // Calculate DV locally (exchange if !exchangeFlux)
      // (p,T), 1/rho, v, e
      if(ComputeDV(threadId))
        return(error::DV);

      if(refRe > 0){

        // Calculate TV locally (exchange if !exchangeFlux)
        // mu, lambda, k ...
        if(ComputeTV(threadId))
          return(error::TV);
      }
      if(refRe > 0 || artificialDissipation){

        if(ComputeVelGrad(threadId))
          return(error::VELGRAD);

        if(artificialDissipation){
          if(ComputeVelDiv(threadId))
            return(error::VELDIV);
        }

        if(ComputeStressTensor(threadId))
          return(error::TAU);
        
        if(ComputeTemperatureGrad(threadId))
          return(error::TGRAD);
        
        if(numScalar > 0){
          if(ComputeConcentrationGrad(threadId)){
            return(error::CONCGRAD);
          }
        }
        // Calculates heat flux (-q_i)
        if(ComputeHeatFlux(threadId))
          return(error::HEATFLUX);
        
        // Computes Q_i = -q_i + (v_j * tau_i_j)
        // stored in heatFluxBuffer
        if(ComputeViscidEnergyFlux(threadId))
          return(error::VEFLUX);
      }


      double m1 = -1.0;
      if(gridType < simulation::grid::RECTILINEAR){
        // Jacobian can be applied in one fell swoop for uniform grids
        std::vector<double> &gridJacobian(gridPtr->Jacobian());
        m1 *= gridJacobian[0];
      }
      
      // For each spatial dimension (iDim):
      for(int iDim = 0;iDim < numDim;iDim++){

        // Populates inviscidFluxBuffer
        if(InviscidFlux(iDim,threadId))
          return(error::IFLUX);

        
        // Populates viscidFluxBuffer
        if(refRe > 0){
          if(ViscousFlux(iDim,threadId))
            return(error::VFLUX);
        }

        if(useWENO){
          // Exchange states
          if(ExchangeState(threadId))
            return(error::XANGE);

          // Exchange fluxes
          if(ExchangeBuffer(fluxMessageId,numEquations,inviscidFluxBuffer,threadId))
            return(error::XANGE);
          
          int wenoStatus = ApplyWENO(iDim,threadId);

          if(wenoStatus > 0){
            // process some WENO error
#pragma omp master
            {
              if(globalPtr){
                std::ostringstream statusMessage;
                statusMessage << "ApplyWENO returned error code (" 
                  << wenoStatus << "), exiting."
                  << std::endl;
                globalPtr->StdOut(statusMessage.str());
                pcpp::io::RenewStream(statusMessage);
              }
            }
            // returns error status to caller
            return(wenoStatus);
          }

          // Sum WENO version of dFluxBuffer to RHS [note use of (m1) uniform gridJac]
          AccumulateToRHS(m1,dFluxBuffer,threadId);

          if(refRe > 0){ // update RHS with viscous part
            if(exchangeFlux){  // exchange flux if required
              if(ExchangeBuffer(fluxMessageId,numEquations,viscidFluxBuffer,threadId))
                return(error::XANGE);
            }
            if(ApplyOperator(iDim,numEquations,viscidFluxBuffer,dFluxBuffer,threadId,true))
              return(error::AO);
            // Sum viscous version of dFluxBuffer to RHS [" note above]
            double m2 = -1.0*m1;
            AccumulateToRHS(m2,dFluxBuffer,threadId,true);
          }
        } else { // no WENO
          if (refRe > 0) { // add viscous flux to inviscid flux
            AccumulateViscousFluxes(-1.0,viscidFluxBuffers,inviscidFluxBuffers,threadId,!exchangeFlux);            
          }
          if(exchangeFlux){  // exchange flux if required
            if(ExchangeBuffer(fluxMessageId,numEquations,inviscidFluxBuffer,threadId))
              return(error::XANGE);
          }

          
          if(ApplyOperator(iDim,numEquations,inviscidFluxBuffer,dFluxBuffer,threadId))
            return(error::AO);
          
          // Sum inviscid+viscous version of dFluxBuffer to RHS [" note above]
          AccumulateToRHS(m1,dFluxBuffer,threadId);

        } // no weno
        
        
      } // Loop over spatial dimensions

      if(gridType > simulation::grid::UNIRECT){ // apply Jacobian if necessary
        GridScaleRHS(threadId);
      }
      
      //  ------- Process other sources ---------
      //   1. boundary sources
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("HandleBCs");
        }
      }      
      int returnCode = HandleBoundaryConditions(threadId);
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("HandleBCs");
        }
      }      
      if(returnCode > 0)
        return(error::BC);
      
      

      //   2. other sources
      returnCode = HandleSources(threadId);


      if(returnCode > 0)
        return(returnCode+200);
      
#pragma omp barrier
      return(0);
    };
    
    // === Secondary interface/wrapper functions =====

    // Computes dFluxBuffer using WENO reconstruction on fluxBuffer
    //
    // Inputs/Outputs
    //  iDim      -- Dimension along which to operate
    //  threadId  -- OpenMP thread number
    int ApplyWENO(int iDim, int threadId) {
#pragma omp barrier
      if(globalPtr){
#pragma omp master 
        {
          globalPtr->FunctionEntry("ApplyWENO");
        }
      }

      int winWidth = coeffsWENO.hiAll - coeffsWENO.loAll + 1;

      double fluxWindow[winWidth*numEquations];
      double projFluxWindow[winWidth*numEquations];
      double projFluxReconst[numEquations];

      double stateL[numEquations];
      double stateR[numEquations];
      double deltaState[numEquations];
      double projDeltaState[numEquations];
      double fluxL[numEquations];
      double fluxR[numEquations];

      double tmat[numEquations*numEquations];
      double tinv[numEquations*numEquations];
      double lambda[numEquations*numEquations];
      double lambdaL[numEquations*numEquations];
      double lambdaR[numEquations*numEquations];
      double normVec[numDim];
      for (int i=0; i<numDim; i++) {
        if (i == iDim) normVec[i] = 1;
        else normVec[i] = 0;
      }

      const std::vector<pcpp::IndexIntervalType> 
        &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());

      const pcpp::IndexIntervalType &regionInterval(threadIntervals[threadId]);

      pcpp::IndexIntervalType bufferInterval;
      bufferInterval.InitSimple(bufferSizes);

      int jDim = (iDim + 1)%numDim;
      int kDim = (iDim + 2)%numDim; // not needed if numDim == 2

      size_t js = regionInterval[jDim].first;
      size_t je = regionInterval[jDim].second;

      // if numDim == 2 we only want the "k" loop to iterate one time
      size_t ks = 0;
      size_t ke = 0;
      if (numDim == 3) {
        ks = regionInterval[kDim].first;
        ke = regionInterval[kDim].second;
      }

      for (size_t k=ks; k<=ke; k++) {
        pcpp::IndexIntervalType opInterval(regionInterval);

        // In the case numDim == 2, kDim == iDim, so we don't want to break
        //  anything by setting opInterval[kDim]
        if (numDim == 3) {
          opInterval[kDim].first = k;
          opInterval[kDim].second = k;
        }

        for (size_t j=js; j<=je; j++) {
          opInterval[jDim].first = j;
          opInterval[jDim].second = j;

          size_t is = regionInterval[iDim].first;
          size_t ie = regionInterval[iDim].second;
          // Note: This routine assumes coeffsWENO.loAll <= 0, coeffsWENO.hiAll >= 0
          opInterval[iDim].first = is + coeffsWENO.loAll - 1;  // -1 because we need to compute WENO flux for left face of leftmost point
          opInterval[iDim].second = ie + coeffsWENO.hiAll;

          std::vector<size_t> indices;
          bufferInterval.GetFlatIndices(opInterval, indices);

          // setup window with fluxes for i == -1
          for (int loc=0; loc<winWidth; loc++) {
            for (int eqn=0; eqn<numEquations; eqn++) {
              //              fluxWindow[eqn*winWidth + loc] = fluxBuffer[eqn*numPointsBuffer + indices[loc]];
              fluxWindow[eqn*winWidth + loc] = inviscidFluxBuffer[eqn*numPointsBuffer + indices[loc]];
            }
          }

          // ii = index into indices, not necessarily falling between is and ie
          for (size_t ii=-coeffsWENO.loAll; ii<indices.size()-coeffsWENO.hiAll; ii++) {
            // get left and right state vectors
            for (int eqn=0; eqn<numEquations; eqn++) {
              stateL[eqn] = myStateBuffers[eqn][indices[ii]];
              stateR[eqn] = myStateBuffers[eqn][indices[ii+1]];
              deltaState[eqn] = stateR[eqn] - stateL[eqn];
            }

            // First use tmat, tinv as temp buffers to get eigenvalues of
            //  Jacobians evaluated at left and right states (needed for
            //  entropy fix computation)
            double gamma = eosPtr->GetGamma();
            FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
              (&numDim, &stateL[0], &stateL[0], &gamma, &normVec[0], &tmat[0], &tinv[0], &lambdaL[0]);
            FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
              (&numDim, &stateR[0], &stateR[0], &gamma, &normVec[0], &tmat[0], &tinv[0], &lambdaR[0]);

            // Now fill tmat, tinv, lambda with Roe eigensystem from
            //  intermediate state
            FC_MODULE(satutil,sat_form_roe_matrices,SATUTIL,SAT_FORM_ROE_MATRICES)
              (&numDim, &stateL[0], &stateR[0], &gamma, &normVec[0], &tmat[0], &tinv[0], &lambda[0]);

            WENO::Project(numEquations, winWidth, &tinv[0], &fluxWindow[0], &projFluxWindow[0]);

            // setup for entropy fix
            WENO::Project(numEquations, 1, &tinv[0], &deltaState[0], &projDeltaState[0]);
            double eta = WENO::EntropyFixEta(numEquations, &lambdaL[0], &lambdaR[0]);

            for (int eqn=0; eqn<numEquations; eqn++) {
              int diag = eqn*(numEquations+1); // get diagonal element of lambda matrix
              int group = (1 + (int) copysign(1.0, lambda[diag]))/2;

              int startLoc = coeffsWENO.LoGroup(group) - coeffsWENO.loAll;
              WENO::ReconstPointVal(&projFluxWindow[eqn*winWidth + startLoc], coeffsWENO, group, projFluxReconst[eqn]);

              double lambdaAbs = std::abs(lambda[diag]);
              double entropyFix = (eta > lambdaAbs ? 0.5*(eta - lambdaAbs) : 0);
              entropyFix *= projDeltaState[eqn];

              //if (entropyFix > 0) {
              //  std::cout << "eta = " << eta << std::endl;
              //  std::cout << "lambdaAbs = " << lambdaAbs << std::endl;
              //  std::cout << "projDeltaState[" << eqn << "] = " << projDeltaState[eqn] << std::endl;
              //  std::cout << "entropyFix[" << eqn << "] = " << entropyFix << std::endl;
              //}

              projFluxReconst[eqn] -= entropyFix;
            }

            WENO::Project(numEquations, 1, &tmat[0], &projFluxReconst[0], &fluxR[0]);

            for (int eqn=0; eqn<numEquations && ii > -coeffsWENO.loAll; eqn++) {
              dFluxBuffer[eqn*numPointsBuffer + indices[ii]] = fluxR[eqn] - fluxL[eqn];
            }

            // shift fluxWindow to right
            for (int loc=1; loc<winWidth; loc++) {
              for (int eqn=0; eqn<numEquations; eqn++) {
                int currIndex = eqn*winWidth + loc;
                int newIndex = eqn*winWidth + (loc-1);
                fluxWindow[newIndex] = fluxWindow[currIndex];
              }
            }

            // copy new values into fluxWindow
            for (int eqn=0; ((eqn<numEquations) && (ii < (indices.size()-coeffsWENO.hiAll-1))); eqn++) {
              //              fluxWindow[eqn*winWidth + (winWidth-1)] = fluxBuffer[eqn*numPointsBuffer + indices[ii+coeffsWENO.hiAll+1]];
              fluxWindow[eqn*winWidth + (winWidth-1)] = inviscidFluxBuffer[eqn*numPointsBuffer + indices[ii+coeffsWENO.hiAll+1]];
            }

            for (int eqn=0; eqn<numEquations; eqn++) {
              fluxL[eqn] = fluxR[eqn];
            }

          } // loop over points in iDim
        } // loop over points in jDim
      } // loop over points in kDim

#pragma omp barrier
      if(globalPtr){
#pragma omp master 
        {
          globalPtr->FunctionExit("ApplyWENO");
        }
      }

      return 0;
    };

    int ArtificialDissipation(int threadId)
    {
      size_t *dOpInterval = &fortranPartitionBufferIntervals[threadId][0];
      std::vector<double>         &gridMetric(gridPtr->Metric());
      std::vector<double>         &gridJacobian(gridPtr->Jacobian());
      
#pragma omp barrier
      if(globalPtr){
#pragma omp master 
        {
          globalPtr->FunctionEntry("ArtificialDissipation");
        }
      }
      
      double dilRamp   = 10.0;
      double sigmaDil  = inputSigmaDil;
      double sigmaDiss = inputSigma;
      double dilCutoff = -1.5;
      
      FC_MODULE(satutil,dissipationweight,SATUTIL,DISSIPATIONWEIGHT)
        (&numDim,&bufferSizes[0],&numPointsBuffer,dOpInterval,&sigmaDiss,
         &sigmaDil,&dilRamp,&dilCutoff,velDiv,sigmaDissipation);

      // Communicate RHS
      for(int iDim = 0;iDim < numDim;iDim++){

        int dissDir = iDim + 1;
        
        // Get Ram's alpha weight
        FC_MODULE(metricops,alphaweight,METRICOPS,ALPHAWEIGHT)
          (&numDim,&numPointsBuffer,&bufferSizes[0],dOpInterval,
           &gridType,&gridMetric[0],&dissDir,alphaDissipation);
        
        // Copy State into fluxBuffer
        for(int iEquation= 0;iEquation < numEquations;iEquation++){
          
          double *iFlux    = &inviscidFluxBuffer[iEquation*numPointsBuffer];
          double *stateBuf = myStateBuffers[iEquation];

          
          FC_MODULE(operators,assignmentyx,OPERATORS,assignmentYX)
            (&numDim,&numPointsBuffer,&bufferSizes[0],dOpInterval,
             stateBuf,iFlux);
        }
        
        ExchangeBuffer(fluxMessageId,numEquations,inviscidFluxBuffer,threadId);

#pragma omp barrier
        
        // Apply dissipation operator on state component:
        ApplyDissipationOperator(iDim,numEquations,inviscidFluxBuffer,dFluxBuffer,false,threadId);
        
        for(int iEquation= 0;iEquation < numEquations;iEquation++){
          
          double *dFlux = &dFluxBuffer[iEquation*numPointsBuffer];
          double *iFlux = &inviscidFluxBuffer[iEquation*numPointsBuffer];
          // Repack the inviscidFluxBuffer by
          // performing R. Vishnampet's metric-based weighting:
          //  ... Malpha=sqrt(M_{alpha,i} * M_{alpha,i}) ; alpha=dir
          //  ... iFlux = alphaDissipation * dFlux
          FC_MODULE(operators,zxy,OPERATORS,ZXY)
            (&numDim,&numPointsBuffer,&bufferSizes[0],dOpInterval,alphaDissipation,
             dFlux,iFlux);
        }

        // Exchange inviscidFluxBuffer for stencil op
        ExchangeBuffer(fluxMessageId,numEquations,inviscidFluxBuffer,threadId);

#pragma omp barrier       

        // Apply split operator on inviscidFluxBuffer:
        ApplyDissipationOperator(iDim,numEquations,inviscidFluxBuffer,dFluxBuffer,true,threadId);

        // Update RHS
        for(int iEquation= 0;iEquation < numEquations;iEquation++){
          
          double *dFlux  = &dFluxBuffer[iEquation*numPointsBuffer];
          double *rhsPtr = myRHSBuffers[iEquation];

          FC_MODULE(operators,ywxpy,OPERATORS,YWXPY)
            (&numDim,&numPointsBuffer,&bufferSizes[0],dOpInterval,
             sigmaDissipation,dFlux,rhsPtr);
          
        }

      } // Loop over spatial dimensions
#pragma omp barrier
      {
        if(globalPtr){
#pragma omp master
          globalPtr->FunctionExit("ArtificialDissipation");
        }
      }
      return(0);
    }; // Artificial dissipation
    
    // No sources (yet)
    int HandleSources(const pcpp::IndexIntervalType &regionInterval)
    {return(0);};
    
    // No sources (yet)
    int HandleSources(int myThreadId){
      int returnCode = 0;
      
      if(artificialDissipation){
        returnCode = ArtificialDissipation(myThreadId);
      }
      
      return(0); 
    };
    
    
    int HandleBoundaryConditions(int threadId)
    {

      // note: profile timers cannot be safely used here
      // due to thread-specific processing performed in this 
      // routine. uncomment bc timers below only for problems 
      // designed to properly measure BC performance.
      
      if(!domainBCsPtr)
        return(0);
      else if (!subRegionsPtr || !domainBoundariesPtr)
        return(1);

      std::vector<BoundaryType>   &domainBoundaries(*domainBoundariesPtr);
      std::vector<BCType>         &domainBCs(*domainBCsPtr);
      std::vector<GridRegionType> &gridRegions(*subRegionsPtr);
      std::vector<double>         &gridMetric(gridPtr->Metric());
      std::vector<double>         &gridJacobian(gridPtr->Jacobian());

      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(6,0.0);

      double gamma = eosPtr->GetGamma();
      gasParams[0] = refRe;
      gasParams[1] = gamma;
      gasParams[2] = specGasConst;
      gasParams[3] = temperatureScale;
      gasParams[4] = refSndSpd;
      gasParams[5] = std::max(gamma/refPr,5.0/3.0);

      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 == navierstokes::HOLE)
            continue;
          
          int normalDir = gridRegion.normalDirection;
          pcpp::IndexIntervalType &regionInterval(gridRegion.regionInterval);
          std::vector<size_t> regionSizes(regionInterval.Sizes());
          size_t numPointsRegion = regionInterval.NNodes();
          pcpp::IndexIntervalType &regionOpBufferInterval(gridRegion.threadPartitionBufferIntervals[threadId]);
          pcpp::IndexIntervalType &regionOpInterval(gridRegion.threadPartitionIntervals[threadId]);

          // bc-associated parameters from input file 
          StateType &bcParamState(boundaryCondition.Params());
          double *bcParams = bcParamState.template GetFieldBuffer<double>("Numbers");
          int *bcFlags     = bcParamState.template GetFieldBuffer<int>("Flags");

          const std::string &bcName(boundaryCondition.BCName());
          const std::string &regionName(gridRegion.regionName);
          const std::string &boundaryName(domainBoundary.Name());
          
          //           std::cout << "bcName " << bcName << " regionName " << regionName << " boundaryName "
          //                     << boundaryName << " bcType " << bcType << std::endl;
          //           int numParams = bcParamState.GetFieldMetaData("Numbers").ncomp;
          //           std::cout << "&bcParams = " << bcParams << std::endl
          //                     << "Number of parameters: " << numParams << std::endl; 
          //             for(int i=0; i<numParams; i++){
          //             std::cout << "i " << i << " bcParams[i] " << bcParams[i] << std::endl;
          //           }
          
          switch(bcType) {
            
          case navierstokes::SPONGE:
            // ** see timer note above
            //            globalPtr->FunctionEntry("Sponge");
            if(normalDir != 0){
              // bcParams[] = {spongeAmplitude,spongePower}
              navierstokes::LinearSponge(bufferSizes,normalDir,regionInterval,regionOpInterval,
                                         regionOpBufferInterval,bcParams,bcFlags,
                                         numEquations,myStateBuffers,targetStateBuffers,myRHSBuffers);
            } else {
              // Sponge blobs are Gaussian
              // bcParams[] = {spongeAmplitude,spongeWidth}
              navierstokes::SpongeZone(bufferSizes,normalDir,regionInterval,regionOpInterval,
                                       regionOpBufferInterval,bcParams,bcFlags,iMask,holeMask,
                                       numEquations,myStateBuffers,targetStateBuffers,myRHSBuffers);
            }
            //            globalPtr->FunctionExit("Sponge");
            break;

          case navierstokes::SAT_FARFIELD:
            //            globalPtr->FunctionEntry("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],
                &gridType,&gridMetric[0],&gridJacobian[0],bcParams,&gasParams[0],
                rhoPtr,rhoVPtr,rhoEPtr,viscidFluxBuffer,&numScalar,scalarPtr,
                rhoRHSPtr,rhoVRHSPtr,rhoERHSPtr,scalarRHSPtr,
                rhoTargetPtr,rhoVTargetPtr,rhoETargetPtr,scalarTargetPtr);

            //            globalPtr->FunctionExit("Farfield");
            break;

          case navierstokes::SAT_NOSLIP_ISOTHERMAL:
            if(refRe > 0){
              if(numScalar > 0) {
                std::ostringstream errMsg;
                errMsg << "EulerRHS::HandleBoundaryConditions: Error SAT_NOSLIP_ISOTHERMAL "
                       << "is not implemented for scalar transport. Aborting." << std::endl;
                if(globalPtr){
                  globalPtr->StdOut(errMsg.str());
                } else {
                  std::cout << errMsg.str();
                }
                return(1);
              }
              bcParams[3]      = boundaryWeight;

              // Call SAT NOSLIP of PlasComCM
              FC_MODULE(satutil,noslip_isothermal,SATUTIL,NOSLIP_ISOTHERMAL)
                ( &numDim,&bufferSizes[0],&numPointsBuffer,&normalDir,&regionSizes[0],
                  &numPointsRegion,&numBoundaryPoints,&threadPartBufferIndices[0],&gridType,
                  &gridMetric[0],&gridJacobian[0],&bcParams[0],&gasParams[0],rhoPtr,rhoVPtr,
                  rhoEPtr,&numScalar,scalarPtr,rhoRHSPtr,rhoVRHSPtr,rhoERHSPtr,
                  scalarRHSPtr,rhoTargetPtr,rhoVTargetPtr,rhoETargetPtr,
                  scalarTargetPtr,muPtr,lambdaPtr);
              
            } else {

              std::ostringstream errMsg;
              errMsg << "EulerRHS::HandleBoundaryConditions:Error: SAT_NOSLIP_ISOTHERMAL "
                     << "not a valid BC for inviscid (Re = 0) flows. Aborting." << std::endl;
              if(globalPtr){
                globalPtr->StdOut(errMsg.str());
              } else {
                std::cout << errMsg.str();
              }
              return(1);
            }

            break;
          case navierstokes::SAT_SLIP_ADIABATIC:
            bcParams[1]      = boundaryWeight;

            FC_MODULE(satutil,slip_adiabatic,SATUTIL,SLIP_ADIABATIC)
              ( &numDim,&bufferSizes[0],&numPointsBuffer,&normalDir,&regionSizes[0],
                &numPointsRegion,&numBoundaryPoints,&threadPartBufferIndices[0],&gridType,
                &gridMetric[0],&gridJacobian[0],bcParams,&gasParams[0],rhoPtr,rhoVPtr,
                rhoEPtr,&numScalar,scalarPtr,rhoRHSPtr,rhoVRHSPtr,rhoERHSPtr,
                scalarRHSPtr,rhoTargetPtr,rhoVTargetPtr,rhoETargetPtr,scalarTargetPtr);
                
            break;
          default:
            std::ostringstream errMsg;
            errMsg << "EulerRHS::HandleBoundaryConditions:Error: Could not resolve BC("
                   << bcName << ") = " << bcType << ". Aborting." << std::endl;
            if(globalPtr){
              globalPtr->StdOut(errMsg.str());
            } else {
              std::cout << errMsg.str() << std::endl;
            }
            return(1);
            break;
            
          }
          
        }
      }
      
      return(0);
      
    };
    
    int ApplyOperator(int dDir,int numComponents,double *uBuffer,double *duBuffer,
                      int threadId,bool skipDensity=false)
    {
          
      size_t dirOffset = dDir*numPointsBuffer; 
      size_t *dOpInterval = &fortranPartitionBufferIntervals[threadId][0];
      
      int one = 1;
      dDir++; // forticate direction

#pragma omp barrier
      if(globalPtr){
#pragma omp master
        {
          globalPtr->FunctionEntry("ApplyOperator");
        }
      }
      
      int startComponent = 0;
      if(skipDensity) startComponent = 1;
      // applyoperator does equation-at-a-time
      for(int iComp = startComponent;iComp < numComponents;iComp++){
        size_t compOffset = iComp*numPointsBuffer;
        FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
          (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
           &dDir,dOpInterval,&numStencils,stencilSizes,
           stencilStarts,&numStencilValues,stencilWeights,
           stencilOffsets,&stencilConn[dirOffset],
           &uBuffer[compOffset],&duBuffer[compOffset]);
      }
      // applyoperatorv does numComponents (i.e. number of equations) 
      // loop on the inside
      // FC_MODULE(operators,applyoperatorv,OPERATORS,APPLYOPERATORV)
      //   (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
      //    &dDir,opInterval,&numStencils,stencilSizes,
      //    stencilStarts,&numStencilValues,stencilWeights,
      //    stencilOffsets,&stencilConn[dirOffset],
      //    &uBuffer[compOffset],&duBuffer[compOffset]);      

#pragma omp barrier
      if(globalPtr){
#pragma omp master
        {         
          globalPtr->FunctionExit("ApplyOperator");
        }
      }

      return(0);
    };
    
    int ApplyDissipationOperator(int dDir,int numComponents,double *uBuffer,double *duBuffer,
                                 bool split = false,int threadId = 0)
    {
          
      size_t dirOffset = dDir*numPointsBuffer; 
      size_t *dOpInterval = &fortranPartitionBufferIntervals[threadId][0];
      
      int one = 1;
      dDir++; // forticate direction

#pragma omp barrier
      if(globalPtr){
#pragma omp master
        {
          globalPtr->FunctionEntry("ApplyDissOperator");
        }
      }
      
      int numStencilsDiss = 0;
      int *stencilSizesDiss = NULL;
      int *stencilStartsDiss = NULL;
      int numStencilValuesDiss = 0;
      double *stencilWeightsDiss = NULL;
      int *stencilOffsetsDiss = NULL;
      int *stencilConnDiss = NULL;

      // Assign dissipation stencil data 
      if(split){
        numStencilsDiss    = artDissOperatorSplit.numStencils;
        stencilSizesDiss   = artDissOperatorSplit.stencilSizes;
        stencilStartsDiss  = artDissOperatorSplit.stencilStarts;
        stencilWeightsDiss = artDissOperatorSplit.stencilWeights;
        stencilOffsetsDiss = artDissOperatorSplit.stencilOffsets;
        stencilConnDiss    = artDissConnSplit;
      } else {
        numStencilsDiss    = artDissOperator.numStencils;
        stencilSizesDiss   = artDissOperator.stencilSizes;
        stencilStartsDiss  = artDissOperator.stencilStarts;
        stencilWeightsDiss = artDissOperator.stencilWeights;
        stencilOffsetsDiss = artDissOperator.stencilOffsets;
        stencilConnDiss    = artDissConn;
      }

      // applyoperator does equation-at-a-time
      for(int iComp = 0;iComp < numComponents;iComp++){

        size_t compOffset = iComp*numPointsBuffer;

        FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
          (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
           &dDir,dOpInterval,&numStencilsDiss,stencilSizesDiss,
           stencilStartsDiss,&numStencilValuesDiss,stencilWeightsDiss,
           stencilOffsetsDiss,&stencilConnDiss[dirOffset],
           &uBuffer[compOffset],&duBuffer[compOffset]);
      }
      // applyoperatorv does numComponents (i.e. number of equations) 
      // loop on the inside
      // FC_MODULE(operators,applyoperatorv,OPERATORS,APPLYOPERATORV)
      //   (&numDim,&bufferSizes[0],&one,&numPointsBuffer,
      //    &dDir,opInterval,&numStencils,stencilSizes,
      //    stencilStarts,&numStencilValues,stencilWeights,
      //    stencilOffsets,&stencilConn[dirOffset],
      //    &uBuffer[compOffset],&duBuffer[compOffset]);      

#pragma omp barrier
      if(globalPtr){
#pragma omp master
        {         
          globalPtr->FunctionExit("ApplyDissOperator");
        }
      }

      return(0);
    };

    void GridScaleRHS(int threadId){
      
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {         
          globalPtr->FunctionEntry("GridScaleRHS");
        }
      }
      
      size_t *dOpInterval = &fortranPartitionBufferIntervals[threadId][0];
      std::vector<double>         &gridJacobian(gridPtr->Jacobian());
      
      //      if(gridType > simulation::grid::UNIRECT){
      for(int iEquation = 0;iEquation < numEquations;iEquation++){
        double *rhsPtr = myRHSBuffers[iEquation];
        FC_MODULE(operators,yxy,OPERATORS,YXY)
          (&numDim,&numPointsBuffer,&bufferSizes[0],dOpInterval,&gridJacobian[0],rhsPtr);
      }
      //       } else {
      //         for(int iEquation = 0;iEquation < numEquations;iEquation++){
      //           double *rhsPtr = myRHSBuffers[iEquation];
      //           FC_MODULE(operators,xax,OPERATORS,XAX)
      //             (&numDim,&numPointsBuffer,&bufferSizes[0],dOpInterval,&gridJacobian[0],rhsPtr);
      //         }
      //       }
      
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {         
          globalPtr->FunctionExit("GridScaleRHS");
        }
      }
    };
    
    
    void AccumulateToRHS(double a,double *dfBuffer,int threadId,bool skipDensity = false) 
    { 
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {         
          globalPtr->FunctionEntry("AccumulateRHS");
        }
      }
        
      size_t *dOpInterval = &fortranPartitionBufferIntervals[threadId][0];
            
      // Y = a*X + Y 
      int startComponent = 0;
      if(skipDensity) startComponent = 1;
      for(int iComp = startComponent;iComp < numEquations;iComp++){
        FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
          (&numDim,&numPointsBuffer,&bufferSizes[0],dOpInterval,&a,
           &dfBuffer[iComp*numPointsBuffer],myRHSBuffers[iComp]);
      }
      
      
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("AccumulateRHS");
        }
      }
      
    };
    

    int ComputeDV(int threadId)
    {

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ComputeDV");
        }
      }
      
      // Need non-forticated intervals here
      const std::vector<pcpp::IndexIntervalType> 
        &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
      
      const pcpp::IndexIntervalType &regionInterval(threadIntervals[threadId]);
      std::vector<size_t> &fortranOpInterval(fortranPartitionBufferIntervals[threadId]);
      
      double zero = 0.0;
      int returnCode = 0;
      // Computed on local partition only 
      // 1/rho, velocity, internal energy
      returnCode = euler::util::ComputeDVBuffer(regionInterval,bufferSizes,
                                                myStateBuffers,dvBuffers);        
        
      // compute pressure and temperature from the eos
      returnCode = eosPtr->ComputePressureBuffer(regionInterval,bufferSizes);        
      returnCode = eosPtr->ComputeTemperatureBuffer(regionInterval,bufferSizes);        

      if(!returnCode && exchangeDV){
        returnCode = ExchangeDV(threadId);
      }

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ComputeDV");
        }
      }

      return(returnCode);

    };

    int ComputeTV(int threadId)
    {
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ComputeTV");
        }
      }
      
      // Need non-forticated intervals here
      const std::vector<pcpp::IndexIntervalType> 
        &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());

      const pcpp::IndexIntervalType &regionInterval(threadIntervals[threadId]);
      std::vector<size_t> &fortranOpInterval(fortranPartitionBufferIntervals[threadId]);
      
      int returnCode = 0;
      // Computed on local partition only 
      if(transportModel == POWER) {
        double heatCapacityCp = eosPtr->GetHeatCapacityCp();
        returnCode = viscid::util::ComputeTVBufferPower(regionInterval,bufferSizes,
                                                  temperaturePtr,tvBuffers,
                                                  beta,power,bulkViscFac,heatCapacityCp,refPr);
      } else {
        std::ostringstream errMsg;
        errMsg << "EulerRHS::RHS:Error: Unknown transport model in EulerRHS. Aborting."
               << std::endl;
        if(globalPtr){
          globalPtr->StdOut(errMsg.str());
          globalPtr->FunctionExit("ComputeTV");
        } else {
          std::cout << errMsg.str();
        }
        
        return(1);
      }

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ComputeTV");
        }
      }
      
      return(returnCode);
    };
    
    int InviscidFlux(int iDim,int threadId)
    {
      // - Inviscid flux - (local if exchangeFlux,+halos otherwise)

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("InviscidFlux");
        }
      }
      
      size_t *fluxOpInterval = NULL;
      if(exchangeFlux){
        fluxOpInterval = &fortranPartitionBufferIntervals[threadId][0];
      } else {
        fluxOpInterval = &fortranBufferIntervals[threadId][0];
      }
      
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ComputeVHat");
        }
      }

      size_t dirOffset = iDim*numPointsBuffer;
      std::vector<double> &gridMetric(gridPtr->Metric());
      double *velocity = &velocityPtr[0];
      int fluxDim = iDim + 1;

      // 0. Get vHat(iDim) = v * Xi_Hat(iDim)
      FC_MODULE(metricops,vhatcomponent,METRICOPS,VHATCOMPONENT)
        (&numDim,&numPointsBuffer,&bufferSizes[0],fluxOpInterval,
         &gridType,&gridMetric[0],&fluxDim,&velocity[0],velHat);

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ComputeVHat");
          globalPtr->FunctionEntry("Flux1D");
        }
      }

      //   1. get f(iDim,1)          = rho*vHat(iDim)
      //   2. get f(iDim,n={2,3,4})  = rho*vHat(iDim)*v(n) + p*delta_ij
      //   3. get f(iDim,numDim+2)   = (rhoE + p)*vHat(iDim)
      FC_MODULE(euler,flux1d,EULER,FLUX1D)
        (&numDim,&numPointsBuffer,&bufferSizes[0],fluxOpInterval,
         &fluxDim,&gridType,&gridMetric[0],rhoPtr,rhoVPtr,rhoEPtr,
         &velHat[0],pressurePtr,inviscidFluxBuffer);

      
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("Flux1D");
        }
      }

      if(numScalar > 0) {

        if(globalPtr){
#pragma omp barrier
#pragma omp master
          {
            globalPtr->FunctionEntry("ScalarFlux1D");
          }
        }
        
        //   4. get f(iDim,numDim+3:+) = rho*Y_n*vHat(iDim) [scalar advection]
        FC_MODULE(euler,scalarflux1d,EULER,SCALARFLUX1D)
          (&numDim,&numPointsBuffer,&bufferSizes[0],fluxOpInterval,
           &numScalar,scalarPtr,&velHat[0],inviscidFluxBuffers[numDim+2]);

        if(globalPtr){
#pragma omp barrier
#pragma omp master
          {
            globalPtr->FunctionExit("ScalarFlux1D");
          }
        }
        
      }

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("InviscidFlux");
        }
      }
      
      return(0);
    };

    int ComputeVelDiv(int threadId) // assumes velgrad is already computed
    {
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ComputeVelDiv");
        }
      }

      const std::vector<pcpp::IndexIntervalType> 
        &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
      const pcpp::IndexIntervalType &regionInterval(threadIntervals[threadId]);
      std::vector<size_t> &fortranOpInterval(fortranPartitionBufferIntervals[threadId]);
      std::vector<double>         &gridMetric(gridPtr->Metric());
      std::vector<double>         &gridJacobian(gridPtr->Jacobian());
      
      // Zero the veldiv buffer [kernel used below is of (+=) type]
      double zero = 0.0;
      FC_MODULE(operators,assignmentxa,OPERATORS,ASSIGNMENTXA)
        (&numDim,&numPointsBuffer,&bufferSizes[0],
         &fortranOpInterval[0],&zero,velDiv);
 
      FC_MODULE(metricops,ijkgradtoxyzdiv,METRICOPS,IJKGRADTOXYZDIV)
        (&numDim,&numPointsBuffer,&bufferSizes[0],&fortranOpInterval[0],
         &gridType,&gridJacobian[0],&gridMetric[0],velGradBufferPtr,
         velDiv);
      
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ComputeVelDiv");
        }
      }

      return(0);

    };

    int ComputeScalarGradient(double *scalarData,double *gradData,int threadId = 0,bool communicateData = false)
    {
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ComputeScalarGrad");
        }
      }
      const size_t *opInterval = &fortranPartitionBufferIntervals[threadId][0];
      std::vector<double>         &gridMetric(gridPtr->Metric());
      std::vector<double>         &gridJacobian(gridPtr->Jacobian());

      if(communicateData){
        ExchangeBuffer(scalarMessageId,1,scalarData,threadId);
      }

      // calculate d(phi)/d(Xi)_i
      for(int i=0;i<numDim;i++){
        size_t gradOffset = numPointsBuffer*i;
        ApplyOperator(i,1,scalarData,&tempVector[gradOffset],threadId);
      }
      
      FC_MODULE(metricops,gradijktogradxyz,METRICOPS,GRADIJKTOGRADXYZ)
        (&numDim,&numPointsBuffer,&bufferSizes[0],opInterval,
         &gridType,&gridJacobian[0],&gridMetric[0],tempVector,gradData);
      
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ComputeScalarGrad");
        }
      }
      
      return(0);
    };
    

    int ComputeVelGrad(int threadId)
    {

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ComputeVelGrad");
        }
      }
      

      // calculate dV_j/d(Xi)_i
      for(int i=0;i<numDim;i++){
        ApplyOperator(i,numDim,velocityPtr,velGradBuffers[i],threadId);
      }

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ComputeVelGrad");
        }
      }
      
      return(0);
    };

    int ComputeConcentrationGrad(int threadId)
    {
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ComputeConcGrad");
        }
      }
      const size_t *opInterval = &fortranPartitionBufferIntervals[threadId][0];

      if(numScalar > 0){
        for(int iScalar = 0;iScalar < numScalar;iScalar++){
          size_t scalarOffset = iScalar*numPointsBuffer;
          size_t scalarGradOffset = numDim*scalarOffset;
          // (1) Y = rhom1 * rhoY
          FC_MODULE(operators,zxy,OPERATORS,ZXY)
            (&numDim,&numPointsBuffer,&bufferSizes[0],opInterval,
             rhom1Ptr,&scalarPtr[scalarOffset],tempScalar);
          // (2) Grad(Y)  = J * g_ij * grad(Y)_i (this is a gradient in physical coordinates)
          ComputeScalarGradient(tempScalar,&scalarGrad[scalarGradOffset],threadId,true);
          //
          // ------------------- NOTES -----------
          // ** This next kernel (YXY) applies the scalar diffusivities
          //    to Grad(Y). Please observe the following:
          // - The scalarDiffusivity should *include* any nondim or scaling factors
          // - ComputeTV or similar should be modified to compute the scalarDiffusivity
          // - Until computed, scalarDiffusivity takes value 1.0 or the user-specified 
          //   value from the input parameter (inputDiffusivity)
          // - BTW: scalarDiffusivity is synomynous with tvBuffers[3+iScalar] 
          // - The scalarDiffusivity is applied *in-place* - which means that the
          //   variable "scalarGrad" now actually contains scalarDiffusivity*Grad(Y).
          // -------------------
          // (3) X_j = scalarDiffusivity*Grad(Y)
          for(int iDim = 0;iDim < numDim; iDim++){
            size_t dirOffset = iDim*numPointsBuffer;
            FC_MODULE(operators,yxy,OPERATORS,YXY)
              (&numDim,&numPointsBuffer,&bufferSizes[0],opInterval,
               &scalarDiffusivity[scalarOffset],&scalarGrad[scalarGradOffset+dirOffset]);
          }
        }
      }

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ComputeConcGrad");
        }
      }
      return(0);
    };
    int ComputeTemperatureGrad(int threadId)
    {
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ComputeTempGrad");
        }
      }
      
      // Need non-forticated intervals here
      const std::vector<pcpp::IndexIntervalType> 
        &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
      const pcpp::IndexIntervalType &regionInterval(threadIntervals[threadId]);
      std::vector<size_t> &fortranOpInterval(fortranPartitionBufferIntervals[threadId]);

      // calculate dT/dX_j
      for(int i=0;i<numDim;i++){
        ApplyOperator(i,1,temperaturePtr,temperatureGradBuffers[i],threadId);
      }
      
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ComputeTempGrad");
        }
      }
      
      return(0);
    };

    int ComputeStressTensor(int threadId)
    {
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ComputeStressTensor");
        }
      }
      
      const std::vector<pcpp::IndexIntervalType> 
        &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
      const pcpp::IndexIntervalType &regionInterval(threadIntervals[threadId]);
      std::vector<double>         &gridMetric(gridPtr->Metric());
      std::vector<double>         &gridJacobian(gridPtr->Jacobian());

      int returnCode = viscid::util::ComputeTauBuffer(regionInterval, bufferSizes, 
                                                      gridType,gridMetric, gridJacobian, refRe, 
                                                      tvBuffers, velGradBuffers, tauBuffers);
      

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ComputeStressTensor");
        }
      }
      
      return(0);
    };

    int ComputeHeatFlux(int threadId)
    {

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ComputeHeatFlux");
        }
      }
      
      const std::vector<pcpp::IndexIntervalType> 
        &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
      const pcpp::IndexIntervalType &regionInterval(threadIntervals[threadId]);
      std::vector<double>         &gridMetric(gridPtr->Metric());
      std::vector<double>         &gridJacobian(gridPtr->Jacobian());

     
      int returnCode = viscid::util::ComputeHeatFluxBuffer(regionInterval, bufferSizes, 
                                                           gridType,gridMetric, gridJacobian, refRe,
                                                           tvBuffers, temperatureGradBuffers, heatFluxBuffers);


      
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ComputeHeatFlux");
        }
      }
      return(0);
    };

    int ComputeViscidEnergyFlux(int threadId)
    {

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ViscidEnergyFlux");
        }
      }
      
      const std::vector<pcpp::IndexIntervalType> 
        &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
      const pcpp::IndexIntervalType &regionInterval(threadIntervals[threadId]);
      const size_t *opInterval = &fortranPartitionBufferIntervals[threadId][0];

      for(int iDim = 0;iDim < numDim;iDim++){
        
        // Compute Q_i = -q_i + (v_j * tau_i_j)
        for(int iVel = 0;iVel < numDim;iVel++){
          int tensorIndex = iDim*numDim + iVel;
          size_t dirOffset = iVel*numPointsBuffer;

          FC_MODULE(operators,ywxpy,OPERATORS,YWXPY)
            (&numDim,&numPointsBuffer,&bufferSizes[0],opInterval,
             &velocityPtr[dirOffset],tauBuffers[tensorIndex],heatFluxBuffers[iDim]);

        }
      }
      
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ViscidEnergyFlux");
        }
      }
      return(0);
    };

    int ViscousFlux(int iDim,int threadId)
    {
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ViscousFlux");
        }
      }
      
      size_t *fluxOpInterval = NULL;
      if(exchangeFlux){
        fluxOpInterval = &fortranPartitionBufferIntervals[threadId][0];
      } else {
        fluxOpInterval = &fortranBufferIntervals[threadId][0];
      }
      std::vector<double>         &gridMetric(gridPtr->Metric());

      // stress tensor indexing
      // 2D, x: 0, 1
      //     y: 1, 2
      // 3D, x: 0, 1, 2
      //     y: 1, 3, 4
      //     z: 2, 4, 5
      int fluxDim = iDim + 1;
      int tau1=iDim*numDim;
      int tau2=tau1+1;
      int tau3=tau2+1;

      
      // viscous fluxes
      FC_MODULE(viscid,strongflux1d,VISCID,STRONGFLUX1D)
        (&numDim,&fluxDim,&bufferSizes[0],&numPointsBuffer,fluxOpInterval,
         &gridType,&gridMetric[0],tauBufferPtr,heatFluxBufferPtr,
         viscidFluxBuffer);
      
      if(numScalar > 0){
        // (4) flux_i = Metric_ij * X_j  (i = iDim)
        for(int iScalar = 0;iScalar < numScalar;iScalar++){
          size_t fluxOffset = numDim+2+iScalar;
          size_t scalarGradOffset = iScalar*numDim*numPointsBuffer;
          FC_MODULE(viscid,scalarflux1d,VISCID,SCALARFLUX1D)
            (&numDim,&fluxDim,&bufferSizes[0],&numPointsBuffer,fluxOpInterval,
             &gridType,&gridMetric[0],&scalarGrad[scalarGradOffset],
             viscidFluxBuffers[fluxOffset]);
        }
      }

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ViscousFlux");
        }
      }
      
      return(0);
    };
    

    void AccumulateViscousFluxes(double a,std::vector<double *> &X,std::vector<double *> &Y,
                                 int threadId,bool fullInterval = false)
    {
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {         
          globalPtr->FunctionEntry("AccumulateVFlux");
        }
      }
      
      size_t *bufferOpInterval = NULL;
      if(fullInterval){
        bufferOpInterval = &fortranPartitionBufferIntervals[threadId][0];
      } else {
        bufferOpInterval = &fortranBufferIntervals[threadId][0];
      }
      
      for(int iComp = 1;iComp < numEquations;iComp++){
        FC_MODULE(operators,yaxpy,OPERATORS,YAXPY)
          (&numDim,&numPointsBuffer,&bufferSizes[0],bufferOpInterval,&a,
           X[iComp],Y[iComp]);
      }
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {         
          globalPtr->FunctionExit("AccumulateVFlux");
        }
      }
      
    }
    
    // ============= Messaging =================
    void PackDVMessage(int threadId)
    {
      if(threadId == 0){
        for(int iDim = 0;iDim < numDV;iDim++)
          haloPtr->PackMessageBuffers(dvMessageId,iDim,1,dvBuffers[iDim]);
      }
      return;
    };

    void PackBufferMessage(int messageId,int numEquations,double *dataBuffer,int threadId)
    {
      if(threadId == 0){
        haloPtr->PackMessageBuffers(messageId,0,numEquations,dataBuffer);
      }
    };
    void UnpackBufferMessage(int messageId,int numEquations,double *dataBuffer,int threadId)
    {
      if(threadId == 0){
        haloPtr->UnPackMessageBuffers(messageId,0,numEquations,dataBuffer);
      }
    }
    void ExchangeBufferMessage(int messageId) // call with threadid == 0 *only* (for now)
    {
      haloPtr->SendMessage(messageId,gridNeighbors,*gridCommPtr);
      haloPtr->ReceiveMessage(messageId,gridNeighbors,*gridCommPtr);
    };
    
    void UnpackDVMessage(int threadId)
    {
      if(threadId == 0){
        for(int iDim = 0;iDim < numDV;iDim++)
          haloPtr->UnPackMessageBuffers(dvMessageId,iDim,1,dvBuffers[iDim]);
      }
      return;
    };

    void ExchangeDVMessage() // call with threadid == 0 *only* (for now)
    {
      haloPtr->SendMessage(dvMessageId,gridNeighbors,*gridCommPtr);
      haloPtr->ReceiveMessage(dvMessageId,gridNeighbors,*gridCommPtr);
    };

    int ExchangeDV(int threadId)
    {
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ExchangeDV");
        }
      }

      PackDVMessage(threadId);
#pragma omp barrier      
#pragma omp master
      {
        ExchangeDVMessage();
      }      
#pragma omp barrier
      UnpackDVMessage(threadId);
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ExchangeDV");
        }
      }
      return(0);
      
    };
    
    int ExchangeBuffer(int messageId,int numEquations,double *dataBuffer,int threadId)
    {
      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionEntry("ExchangeBuffer");
        }
      }
      PackBufferMessage(messageId,numEquations,dataBuffer,threadId);
#pragma omp barrier      
      
#pragma omp master
      {
        ExchangeBufferMessage(messageId);
      }
      
#pragma omp barrier
      
      UnpackBufferMessage(messageId,numEquations,dataBuffer,threadId);

      if(globalPtr){
#pragma omp barrier
#pragma omp master
        {
          globalPtr->FunctionExit("ExchangeBuffer");
        }
      }

      return(0);

    };
    
    
    int ExchangeInviscidFlux(int threadId)
    {
      return(ExchangeBuffer(fluxMessageId,numEquations,inviscidFluxBuffer,threadId));
    };

    int ExchangeState(int threadId)
    {
      if(globalPtr){
 #pragma omp barrier
 #pragma omp master
         {
           globalPtr->FunctionEntry("ExchangeState");
         }
      }

#pragma omp master
      {
        haloPtr->PackMessageBuffers(fluxMessageId,0,1,rhoPtr);
        haloPtr->PackMessageBuffers(fluxMessageId,1,numDim,rhoVPtr);
        haloPtr->PackMessageBuffers(fluxMessageId,numDim+1,1,rhoEPtr);

        haloPtr->SendMessage(fluxMessageId,gridNeighbors,gridPtr->Communicator());
        haloPtr->ReceiveMessage(fluxMessageId,gridNeighbors,gridPtr->Communicator());

        haloPtr->UnPackMessageBuffers(fluxMessageId,0,1,rhoPtr);
        haloPtr->UnPackMessageBuffers(fluxMessageId,1,numDim,rhoVPtr);
        haloPtr->UnPackMessageBuffers(fluxMessageId,numDim+1,1,rhoEPtr);
      }

       if(globalPtr){
 #pragma omp barrier
 #pragma omp master
         {
           globalPtr->FunctionExit("ExchangeState");
         }
       }

      return 0;
    };

    // ===== State-related interfaces ====

    // *NOT* thread safe - master only
    // 
    int SetParamState(const StateType &paramState)
    {

#pragma omp barrier
      if(globalPtr){
#pragma omp master
        globalPtr->FunctionEntry("SetParamState");
      }
      
      //
      // get non-dimensional reference parameters
      //

      {
        int paramHandle = paramState.GetDataIndex("refRe");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType &paramData(paramState.Field(paramHandle));
          const double *paramPtr = paramData.Data<double>();
          refRe = *paramPtr;
          if(refRe < 0)
            refRe = 0;
        } else {
          refRe = 0.0;
        }
      }

      {
        int paramHandle = paramState.GetDataIndex("inputDiffusivity");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType &paramData(paramState.Field(paramHandle));
          inputDiffusivity = paramData.Data<double>();
        } else {
          inputDiffusivity = NULL;
        }
      }

      {
        int paramHandle = paramState.GetDataIndex("refPr");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType &paramData(paramState.Field(paramHandle));
          const double *paramPtr = paramData.Data<double>();
          refPr = *paramPtr;
          if(refPr < 0)
            refPr = 0.75;
        } else {
          refPr = 0.75;
        }
      }

      //
      // Define the reference state based on the density, pressure, and temperature
      // all other reference parameters derive as follows for a calorically perfect, ideal gas
      //
      // p_ref = p_inf
      // rho_ref = rho_inf
      // temperature_ref = temperature_inf
      //
      // The user has the option to run in dimensional or non-dimensional mode
      // dimensional mode calculates required EOS parameters (i.e. the specific gas constant
      // for a perfect gas) and sets the above to one
      //
      // The default behavior to run in non-dimensional mode with STP reference conditions for air
      //
      {
        int paramHandle = paramState.GetDataIndex("nonDimensional");
        std::ostringstream statusMsg;
        statusMsg << (paramHandle >= 0 ? "User chose" : "Defaulting to"); 
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType &paramData(paramState.Field(paramHandle));
          const int *paramPtr = paramData.Data<int>();
          nonDimensionalMode = *paramPtr;
          statusMsg << (nonDimensionalMode==1 ? " non-dimensional " : 
                        nonDimensionalMode==2 ? " PlasComCM-non-dimensional " :
                        " dimensional ");
        } else {
          nonDimensionalMode = 1;
          statusMsg << " non-dimensional ";
        }
        statusMsg << "mode." << std::endl;
        if(globalPtr){
          globalPtr->StdOut(statusMsg.str());
        } else {
          std::cout << statusMsg.str();
        }
      }

      //
      // get the reference state from the parameters
      // if it is not defined, or incomplete, we specify a default state
      // and warn the user
      //
      bool haveRefRho = false;
      bool haveRefPressure = false;
      bool haveRefTemperature = false;
      {
        int paramHandle = paramState.GetDataIndex("refRho");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType 
            &paramData(paramState.Field(paramHandle));
          const double *paramPtr = paramData.Data<double>();
          refRho = *paramPtr;
          haveRefRho = true;
        }
      }

      {
        int paramHandle = paramState.GetDataIndex("refPressure");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType 
            &paramData(paramState.Field(paramHandle));
          const double *paramPtr = paramData.Data<double>();
          refPressure = *paramPtr;
          haveRefPressure = true;
        }
      }

      {
        int paramHandle = paramState.GetDataIndex("refTemperature");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType 
            &paramData(paramState.Field(paramHandle));
          const double *paramPtr = paramData.Data<double>();
          refTemperature = *paramPtr;
          haveRefTemperature = true;
        }
      }

      if(!(haveRefRho && haveRefPressure && haveRefTemperature)){
        refRho = 1.225;
        refPressure = 101325;
        refTemperature = 288.15;
        std::ostringstream statusMessage;
        statusMessage << "Warning: EOS reference state missing or incomplete in input. " << std::endl
                      << "User must specify refRho, refPressure, and refTemperature. " << std::endl
                      << "Setting to default values: "  << std::endl
                      << "\t refRho=" << refRho << std::endl
                      << "\t refPressure=" << refPressure << std::endl
                      << "\t refTemperature=" << refTemperature << std::endl
                      << std::endl;
        if(globalPtr){
          globalPtr->StdOut(statusMessage.str());
        } else {
          std::cout << statusMessage.str();
        }
      }

      {
        int paramHandle = paramState.GetDataIndex("refLength");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType 
            &paramData(paramState.Field(paramHandle));
          const double *paramPtr = paramData.Data<double>();
          refLength = *paramPtr;
        } else {
          refLength = 1.0;
          std::ostringstream statusMessage;
          statusMessage << "Warning: Could not find refLength in input" << std::endl
                        << "Setting to default value " << refLength << std::endl;
          if(globalPtr){
            globalPtr->StdOut(statusMessage.str());
          } else {
            std::cout << statusMessage.str();
          }
        }
        
        
      }
      
      //
      // initialize equation of state object
      // get data associated with the chosen equation of state (i.e. gamma)
      // and transport laws that may be needed to fully define the reference state
      //
      if(gasModel == IDEAL) {
        eosPtr = new eos::perfect_gas();

        gammaHandle = paramState.GetDataIndex("gamma");
        if(gammaHandle < 0) {
          std::ostringstream statusMessage;
          statusMessage << "Gamma does not exist in input parameter state." 
                        << " Exiting." << std::endl;
          if(globalPtr){
            globalPtr->StdOut(statusMessage.str());
            globalPtr->FunctionExit("SetParamState");
          } else {
            std::cout << statusMessage.str();
          }
          return(1);
        }

        const pcpp::field::dataset::DataBufferType &gammaData(paramState.Field(gammaHandle));
        const double *gammaPtr = gammaData.Data<double>();
        if(!gammaPtr) {
          std::ostringstream statusMessage;
          statusMessage << "Gamma not defined in input parameter state." 
                        << " Exiting." << std::endl;
          if(globalPtr){
            globalPtr->StdOut(statusMessage.str());
            globalPtr->FunctionExit("SetParamState");
          } else {
            std::cout << statusMessage.str();
          }
          return(1);
        }

        eosPtr->SetGamma(*gammaPtr);
      } else {
        std::ostringstream statusMessage;
        statusMessage << "Error: Unknown gas model selected." << gasModel 
                      << " Exiting." << std::endl;
        if(globalPtr){
          globalPtr->StdOut(statusMessage.str());
          globalPtr->FunctionExit("SetParamState");
        } else {
          std::cout << statusMessage.str();
        }
        return(1);
      }

      if(refRe > 0) {
        if(transportModel == POWER) {

          // power
          powerHandle = paramState.GetDataIndex("power");
          if(powerHandle < 0) {
            std::ostringstream statusMessage;
            statusMessage << "Power does not exist in input parameter state." << std::endl
                          << " Required for the Power Law transport model." << std::endl
                          << " Exiting." << std::endl;
            if(globalPtr){
              globalPtr->StdOut(statusMessage.str());
              globalPtr->FunctionExit("SetParamState");
            } else {
              std::cout << statusMessage.str();
            }
            return(1);
          }

          const pcpp::field::dataset::DataBufferType 
            &powerData(paramState.Field(powerHandle));
          powerPtr = powerData.Data<double>();
          if(!powerPtr) {
            std::ostringstream statusMessage;
            statusMessage << "Power not defined in input parameter state."  << std::endl
                          << "Required for the Power Law transport model." << std::endl
                          << "Example: power=2/3 for air." << std::endl
                          << "Exiting." << std::endl;
            if(globalPtr){
              globalPtr->StdOut(statusMessage.str());
              globalPtr->FunctionExit("SetParamState");
            } else {
              std::cout << statusMessage.str();
            }
            return(1);
          }
          power = *powerPtr;

          // beta
          betaHandle = paramState.GetDataIndex("beta");
          if(betaHandle < 0) {
            std::ostringstream statusMessage;
            statusMessage << "beta does not exist in input parameter state." << std::endl
                          << "Required for the beta Law transport model." << std::endl
                          << "Exiting." << std::endl;
            if(globalPtr){
              globalPtr->StdOut(statusMessage.str());
              globalPtr->FunctionExit("SetParamState");
            } else {
              std::cout << statusMessage.str();
            }
            return(1);
          }

          const pcpp::field::dataset::DataBufferType 
            &betaData(paramState.Field(betaHandle));
          betaPtr = betaData.Data<double>();
          if(!betaPtr) {
            std::ostringstream statusMessage;
            statusMessage << "beta not defined in input parameter state." << std::endl
                          << "Required for the beta Law transport model." << std::endl
                          << "Example: beta=4.093*10^-7 for air." << std::endl
                          << "Exiting." << std::endl;
            if(globalPtr){
              globalPtr->StdOut(statusMessage.str());
              globalPtr->FunctionExit("SetParamState");
            } else {
              std::cout << statusMessage.str();
            }
            return(1);
          }
          beta = *betaPtr;

          //bulkViscosity Factor
          bulkViscFacHandle = paramState.GetDataIndex("bulkViscFac");
          if(bulkViscFacHandle < 0) {
            std::ostringstream statusMessage;
            statusMessage << "bulkViscFac does not exist in input parameter state." << std::endl
                          << "Required for the bulkViscFac Law transport model."<< std::endl
                          << "Exiting." << std::endl;
            if(globalPtr){
              globalPtr->StdOut(statusMessage.str());
              globalPtr->FunctionExit("SetParamState");
            } else {
              std::cout << statusMessage.str();
            }
            return(1);
          }

          const pcpp::field::dataset::DataBufferType 
            &bulkViscFacData(paramState.Field(bulkViscFacHandle));
          bulkViscFacPtr = bulkViscFacData.Data<double>();
          if(!bulkViscFacPtr) {
            std::ostringstream statusMessage;
            statusMessage << "bulkViscFac not defined in input parameter state." << std::endl
                          << "Required for the bulkViscFac Law transport model."<< std::endl
                          << "Example: bulkViscFac=0.6 for air." << std::endl
                          << "Exiting." << std::endl;
            if(globalPtr){
              globalPtr->StdOut(statusMessage.str());
              globalPtr->FunctionExit("SetParamState");
            } else {
              std::cout << statusMessage.str();
            }
            return(1);
          }
          bulkViscFac = *bulkViscFacPtr;
          
        } else {
          std::ostringstream statusMessage;
          statusMessage << "Error: Unknown transport model selected." 
                        << transportModel << " Exiting." << std::endl;
          if(globalPtr){
            globalPtr->StdOut(statusMessage.str());
            globalPtr->FunctionExit("SetParamState");
          } else {
            std::cout << statusMessage.str();
          }
          return(1);
        }
      }
      
      //
      // derive the remaining reference state variables
      //
      if(gasModel == IDEAL) {
        double gamma = eosPtr->GetGamma();
        if(nonDimensionalMode == 1){
          specGasConst = 1.0;
          refSpecGasConst = refPressure/refRho/refTemperature;
          refVelocity = sqrt(refPressure/refRho);
          refSndSpd = refVelocity*sqrt(gamma);
          refEnergy = refPressure/(gamma-1.0)/refRho;
          refCp = refSpecGasConst*gamma/(gamma-1.0);
          refCv = refCp-refSpecGasConst;


          pressureScale = refPressure;
          temperatureScale = refTemperature;
          velocityScale = refVelocity;

        } else if(nonDimensionalMode == 2){

          specGasConst = (gamma-1)/gamma;
          refSpecGasConst = refPressure/refRho/refTemperature;
          refVelocity = sqrt(gamma*refPressure/refRho);
          refSndSpd = refVelocity;
          refEnergy = refPressure/(gamma-1.0)/refRho;
          refCp = refSpecGasConst*gamma/(gamma-1.0);
          refCv = refCp-refSpecGasConst;

          pressureScale = refRho*refSndSpd*refSndSpd;
          temperatureScale = (gamma-1)*refTemperature;
          velocityScale = refSndSpd;

          std::ostringstream statusMessage;
          statusMessage << "Enabling legacy PlasComCM nonDimensionalization mode." << std::endl;
          if(globalPtr){
            globalPtr->StdOut(statusMessage.str());
          } else {
            std::cout << statusMessage.str();
          }

        } else {

          refVelocity = sqrt(refPressure/refRho);
          specGasConst = refPressure/refRho/refTemperature;
          refSndSpd = refVelocity*sqrt(gamma);
          refCp = specGasConst*gamma/(gamma-1.0);
          refCv = refCp-specGasConst;
          refEnergy = refPressure/(gamma-1.0)/refRho;

          refRho = 1.0;
          refPressure = 1.0;
          refTemperature = 1.0;
          refVelocity = 1.0;
          refSpecGasConst = 1.0;

          pressureScale = 1.0;
          temperatureScale = 1.0;
          velocityScale = 1.0;
        }

        eosPtr->SetSpecificGasConstant(specGasConst);
        eosPtr->InitializeMaterialProperties();
      }


      //
      // in dimensional mode all the non-dimensional numbers become unity
      // 
      if(refRe>0){
        if(nonDimensionalMode != 0){
          refMu = refRho*refVelocity*refLength/refRe;
          refKappa = refCp*refMu/refPr;
        } else {
          refRe = 1.0;
          refMu = 1.0;
        }
      }

      //
      // Scale the input parameters for non-dimensional runs
      //
      if(refRe>0){
        if(nonDimensionalMode == 1){
          if(transportModel == POWER){
            beta = 1.0;
          }
        } else if(nonDimensionalMode == 2){
          double gamma = eosPtr->GetGamma();
          if(transportModel == POWER){
            //beta = beta*pow(refTemperature,power)/refMu;
            beta = pow(gamma-1,power);
          }
        } 
      }


      //
      // Report the current reference state to user
      std::ostringstream refStateMsg;
      refStateMsg << "EulerRHS Reference State: "          << std::endl
                  << "Reynolds Number: \t"        << refRe << std::endl
                  << "Prandtl Number: \t"         << refPr << std::endl
        
                  << "Length Reference: \t"       << refLength      << std::endl
                  << "Density Reference: \t"      << refRho         << std::endl
                  << "Pressure Reference: \t"     << refPressure    << std::endl
                  << "Temperature Reference: \t"  << refTemperature << std::endl

                  << "Specific Gas Constant Reference: \t" << refSpecGasConst << std::endl
                  << "Energy Reference: \t"                << refEnergy       << std::endl
                  << "Velocity Reference: \t"              << refVelocity     << std::endl
                  << "Speed of Sound Reference: \t"        << refSndSpd       << std::endl
                  << "Cp Reference: \t"                    << refCp           << std::endl
                  << "Cv Reference: \t"                    << refCv           << std::endl;
        
      if(nonDimensionalMode == 0) {
        refStateMsg << "EulerRHS set to run in dimensional mode."    << std::endl
                    << "Specific Gas Constant: \t" << specGasConst   << std::endl
                    << "Speed of sound: \t"        << refSndSpd      << std::endl
                    << "Cp: \t"                    << eosPtr->GetHeatCapacityCp() << std::endl
                    << "Cv: \t"                    << eosPtr->GetHeatCapacityCv() << std::endl;
      }
      
      if(globalPtr){
        globalPtr->StdOut(refStateMsg.str());
      } else {
        std::cout << refStateMsg.str();
      }
      
      std::ostringstream paramStatusMsg;

      //
      // initialization of other simulation parameters
      //
      inputCFL = 1.0;
      int inputCFLHandle = paramState.GetDataIndex("inputCFL");
      if(inputCFLHandle >= 0){
        const pcpp::field::dataset::DataBufferType &cflData(paramState.Field(inputCFLHandle));
        const double *inCFLPtr = cflData.Data<double>();
        inputCFL = *inCFLPtr;
        paramStatusMsg << "Input CFL = " << inputCFL << std::endl;
      }
      if(inputCFL <= 0){
        inputCFL = 1.0;
        paramStatusMsg << "Defaulting CFL = " << inputCFL << std::endl;
      }
      
      
      inputDT = -1.0;
      int inputDTHandle = paramState.GetDataIndex("inputDT");
      if(inputDTHandle  >= 0){
        const pcpp::field::dataset::DataBufferType &inputDTData(paramState.Field(inputDTHandle));
        const double *inputDTPtr = inputDTData.Data<double>();
        inputDT = *inputDTPtr;
        paramStatusMsg << "Input DT = " << inputDT << std::endl;
      } 
      if(inputDT < 0){ // inputDT = 0, rhs and advancer still operate, but state doesn't change (useful for testing)
        timeSteppingMode = CONSTANT_CFL;
        paramStatusMsg << "User-specified DT is invalid, defaulting to CONSTANT_CFL mode."
                       << std::endl; 
      }
      
      sigmaHandle = paramState.GetDataIndex("sigmaDissipation");
      if(sigmaHandle >= 0){
        const pcpp::field::dataset::DataBufferType &sigmaData(paramState.Field(sigmaHandle));
        sigmaPtr = sigmaData.Data<double>();
        inputSigma = *sigmaPtr;
      } else {
        sigmaPtr = NULL;
        inputSigma = 0.0;
      }
      
      int sigmaDilHandle = paramState.GetDataIndex("sigmaDilatation");
      if(sigmaDilHandle >= 0){
        const pcpp::field::dataset::DataBufferType &sigmaData(paramState.Field(sigmaDilHandle));
        sigmaDilPtr = sigmaData.Data<double>();
        inputSigmaDil = *sigmaDilPtr;
      } else {
        sigmaDilPtr = NULL;
        inputSigmaDil = 0.0;
      }


      if(inputSigma > 0.0 || inputSigmaDil > 0.0){
        artificialDissipation = true;
      }
      if(artificialDissipation){
        if(inputSigma < 0 || inputSigmaDil < 0){
          return(1);
        }
      }
      


      {
        int paramHandle = paramState.GetDataIndex("Flags");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType &paramData(paramState.Field(paramHandle));
          const int *paramPtr = paramData.Data<int>();
          int optionFlags = *paramPtr;
          if(*paramPtr > 0){
            //            useAlternateRHS = optionFlags&ALTERNATERHS;
            //            exchangeFlux    = optionFlags&EXCHANGEFLUX;
            exchangeDV      = optionFlags&EXCHANGEDV;
            useWENO         = optionFlags&USEWENO;
          }
        } 
        if(refRe > 0){
          paramStatusMsg << "EulerRHS: Viscous flow, enabling DV exchange." << std::endl;
          exchangeDV = true;
        }
        if(artificialDissipation){
          paramStatusMsg << "EulerRHS: Artificial dissipation ON.";
          if(refRe == 0)
            paramStatusMsg << " Enabling DV exchange.";
          paramStatusMsg << std::endl;
          exchangeDV = true;
        }
        if(exchangeFlux){
          paramStatusMsg << "EulerRHS: Enabling flux exchange." << std::endl;
        }
      }

#pragma omp barrier
      if(globalPtr){
#pragma omp master
        {
          globalPtr->StdOut(paramStatusMsg.str());
          globalPtr->FunctionExit("SetParamState");
        }
      } else {
        std::cout << paramStatusMsg.str();
      }

      return(0);
    };
    
    // *NOT* thread-safe, call on master only
    int InitializeState(StateType &inState)
    {      

      if(!gridInitialized)
        return(1);
      
      baseStatePtr = &inState;

      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;

      cflHandle = inState.GetDataIndex("simCFL");
      if(cflHandle >= 0){
        pcpp::field::dataset::DataBufferType &cflData(inState.Field(cflHandle));
        cflPtr = cflData.Data<double>();
      }

      dtHandle = inState.GetDataIndex("simDT");
      if(dtHandle >= 0){
        pcpp::field::dataset::DataBufferType &dtData(inState.Field(dtHandle));
        dtPtr = dtData.Data<double>();
      }
      
      //      rhoHandle = inState.GetDataIndex("rho");
      rhoHandle = inState.GetFieldHandle("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.GetFieldDataByHandle(rhoHandle));
      rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      //      rhoVHandle = inState.GetDataIndex("rhoV");
      rhoVHandle = inState.GetFieldHandle("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.GetFieldDataByHandle(rhoVHandle));
      rhoVPtr   = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
        
      //      rhoEHandle = inState.GetDataIndex("rhoE");
      rhoEHandle = inState.GetFieldHandle("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.GetFieldDataByHandle(rhoEHandle));
      rhoEPtr   = rhoEData.Data<double>();
      if(!rhoEPtr)
        return(1);
        
      pressureHandle = inState.GetDataIndex("pressure");
      if(pressureHandle >= 0){
        pcpp::field::dataset::DataBufferType &pressureData(inState.Field(pressureHandle));
        pressurePtr = pressureData.Data<double>();
      }
        
      temperatureHandle = inState.GetDataIndex("temperature");
      if(temperatureHandle >= 0){
        pcpp::field::dataset::DataBufferType &temperatureData(inState.Field(temperatureHandle));
        temperaturePtr = temperatureData.Data<double>();
      }
        
      // Optional data
      rhom1Handle = inState.GetDataIndex("rhom1");
      if(rhom1Handle >= 0){
        pcpp::field::dataset::DataBufferType &rhom1Data(inState.Field(rhom1Handle));
        rhom1Ptr = rhom1Data.Data<double>();
      }

      velocityHandle = inState.GetDataIndex("velocity");
      if(velocityHandle >= 0){
        pcpp::field::dataset::DataBufferType &velocityData(inState.Field(velocityHandle));
        velocityPtr = velocityData.Data<double>();
      }

      internalEnergyHandle = inState.GetDataIndex("internalEnergy");
      if(internalEnergyHandle >= 0){
        pcpp::field::dataset::DataBufferType &internalEnergyData(inState.Field(internalEnergyHandle));
        internalEnergyPtr = internalEnergyData.Data<double>();
      }

      //      scalarHandle = inState.GetDataIndex("scalarVars");
      scalarHandle = inState.GetFieldHandle("scalarVars");
      if(scalarHandle >= 0){
        int scalarDataIndex = inState.GetFieldIndexByHandle(scalarHandle);
        pcpp::field::dataset::DataBufferType &scalarData(inState.GetFieldDataByHandle(scalarHandle));
        numScalar = inState.Meta()[scalarDataIndex].ncomp;
        scalarPtr   = scalarData.Data<double>();
        if(!scalarPtr)
          return(1);
      } else {
        scalarPtr = NULL;
        numScalar = 0;
      }


      int scalarGradHandle = inState.GetDataIndex("scalarGrad");
      if(scalarGradHandle >= 0){
        pcpp::field::dataset::DataBufferType &gradData(inState.GetFieldDataByHandle(scalarGradHandle));
        scalarGrad   = gradData.Data<double>();
        if(!scalarGrad)
          return(1);
      } else {
        scalarGrad = NULL;
      }
      
      int stencilConnHandle = inState.GetDataIndex("sconn");
      if(stencilConnHandle >= 0){
        //         if(!stencilConn) 
        //           return(1);
        //         inState.SetFieldBuffer("sconn",stencilConn);
        pcpp::field::dataset::DataBufferType &scData(inState.Field(stencilConnHandle));
        stencilConn   = scData.Data<int>();
        if(!stencilConn)
          return(1);
      } else {
        stencilConn = NULL;
      }        
      
      int iMaskHandle = inState.GetDataIndex("imask");
      if(iMaskHandle >= 0){
        pcpp::field::dataset::DataBufferType &imData(inState.Field(iMaskHandle));
        iMask   = imData.Data<int>();
        if(!iMask)
          return(1);
      } else {
        iMask = NULL;
      }        

      int dilatationHandle = inState.GetDataIndex("divV");
      if(dilatationHandle >= 0){
        pcpp::field::dataset::DataBufferType &dilData(inState.Field(dilatationHandle));
        velDiv = dilData.Data<double>();
        if(!velDiv)
          return(1);
      } else {
        velDiv = NULL;
      }

      // dynamic viscosity, mu
      muHandle = inState.GetDataIndex("mu");
      if(muHandle >= 0){
        pcpp::field::dataset::DataBufferType &muData(inState.Field(muHandle));
        muPtr = muData.Data<double>();
      } else {
        muPtr = NULL;
      }

      // bulk viscosity, lambda
      lambdaHandle = inState.GetDataIndex("lambda");
      if(lambdaHandle >= 0 ){
        pcpp::field::dataset::DataBufferType &lambdaData(inState.Field(lambdaHandle));
        lambdaPtr = lambdaData.Data<double>();
      } else {
        lambdaPtr = NULL;
      }

      // heat transfer coefficient, kappa
      kappaHandle = inState.GetDataIndex("kappa");
      if(kappaHandle >= 0){
        pcpp::field::dataset::DataBufferType &kappaData(inState.Field(kappaHandle));
        kappaPtr = kappaData.Data<double>();
      } else {
        kappaPtr = NULL;
      }
        
      int diffusivityHandle = inState.GetDataIndex("scalarDiffusivity");
      if(diffusivityHandle >= 0){
        pcpp::field::dataset::DataBufferType &diffusivityData(inState.Field(diffusivityHandle));
        int numDiff = inState.Meta()[diffusivityHandle].ncomp;
        if(numDiff != numScalar){
          std::cerr << "Scalar Diffusivities of improper size ( "
                    << numDiff << " != number of scalars ("
                    << numScalar << ")." << std::endl;
          return(1);
        }
        scalarDiffusivity = diffusivityData.Data<double>();
      } else {
        scalarDiffusivity = NULL;
      }

      // Number of equations we actually solve
      numEquations = numDim+2+numScalar;
      myStateBuffers.resize(numEquations+1,NULL); // extra one for scalar transport to keep API happy
        
      //        myStateBuffers.resize(4,NULL);
      int iEquation = 0;
      myStateBuffers[iEquation++] = rhoPtr;
      for(int iDim = 0;iDim < numDim;iDim++)
        myStateBuffers[iEquation++] = &rhoVPtr[iDim*numPointsBuffer];
      myStateBuffers[iEquation++] = rhoEPtr;
      for(int iScalar = 0;iScalar < numScalar;iScalar++)
        myStateBuffers[iEquation++] = &scalarPtr[iScalar*numPointsBuffer];
      myRHSBuffers.resize(numEquations+1,NULL);

      if(dvBuffers.empty()){

        // 
        numDV = numDim+4; // p T 1/rho V e

        // New shiny
        // dvBuffer = new double [numPointsBuffer * numDV];
        // fluxComponentPtr = new double [numFluxComponents * numPointsBuffer];
        
        // Old busted
        dvBuffers.resize(numDV);   //  = new double [6*numPoints];
        
        //        scaledPressure = new double [numPointsBuffer];

        if(pressurePtr){
          dvBuffers[0] = pressurePtr;
        } else {
          dvBuffers[0] = new double [numPointsBuffer];
          pressurePtr = dvBuffers[0];
        }
        if(temperaturePtr){
          dvBuffers[1] = temperaturePtr;
        } else {
          dvBuffers[1] = new double [numPointsBuffer];
          temperaturePtr = dvBuffers[1];
        }
        if(rhom1Ptr){
          dvBuffers[2] = rhom1Ptr;
        } else {
          dvBuffers[2] = new double [numPointsBuffer];
          rhom1Ptr = dvBuffers[2];
        }
        if(!velocityPtr)
          velocityPtr = new double [numDim*numPointsBuffer];
        for(int iDim = 0;iDim < numDim;iDim++)
          dvBuffers[iDim+3] = &velocityPtr[iDim*numPointsBuffer];
        if(internalEnergyPtr){
          dvBuffers[3+numDim] = internalEnergyPtr;
        } else {
          dvBuffers[3+numDim] = new double [numPointsBuffer];
          internalEnergyPtr = dvBuffers[3+numDim];
        }
      } 

      if(tvBuffers.empty()){

        numTV = 3+2*numScalar; // mu lambda kappa [scalar diffusivities]
        tvBuffers.resize(numTV,NULL);
        scaledTau = new double [numPointsBuffer];

        if(muPtr){
          tvBuffers[0] = muPtr;
        } else {
          tvBuffers[0] = new double [numPointsBuffer];
          muPtr = tvBuffers[0];
        }
        if(lambdaPtr){
          tvBuffers[1] = lambdaPtr;
        } else {
          tvBuffers[1] = new double [numPointsBuffer];
          lambdaPtr = tvBuffers[1];
        }
        if(kappaPtr){
          tvBuffers[2] = kappaPtr;
        } else {
          tvBuffers[2] = new double [numPointsBuffer];
          kappaPtr = tvBuffers[2];
        }
        if(numScalar > 0){
          if(!scalarDiffusivity){
            scalarDiffusivity = new double [numScalar*numPointsBuffer];
          }
          for(int iScalar = 0;iScalar < numScalar;iScalar++){
            size_t scalarOffset = iScalar*numPointsBuffer;
            tvBuffers[3+iScalar] = &scalarDiffusivity[scalarOffset];
            double defaultDiffusivity = 1.0;
            if(inputDiffusivity){
              defaultDiffusivity = inputDiffusivity[iScalar];
            } 
            for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
              scalarDiffusivity[scalarOffset+iPoint] = defaultDiffusivity;
          }

          // tvBuffers[3+numScalar+iScalar] if needed, allocate in the above loop

        }
      }

      // initialize the RHS state object
      rhsState.AddField("vHat",'n',1,8,""); // only one component of vhat
      rhsState.AddField("viscidEnergy",'n',3,8,"");
      rhsState.AddField("inviscidFlux",'n',numEquations,8,"");
      rhsState.AddField("viscidFlux",'n',numEquations,8,"");
      rhsState.AddField("dFlux",'n',numEquations,8,"");
      if(refRe > 0 || artificialDissipation){
        rhsState.AddField("velGrad",'n',numDim*numDim,8,"");
        if(artificialDissipation){
          rhsState.AddField("velDiv",'n',1,8,"");
          rhsState.AddField("sigmaDiss",'n',1,8,"");
          rhsState.AddField("alphaDiss",'n',1,8,"");
        }
      }
      if(refRe > 0){
        rhsState.AddField("tau",'n',(numDim*(numDim+1))/2,8,"");
        rhsState.AddField("temperatureGrad",'n',numDim,8,"");
        rhsState.AddField("heatFlux",'n',numDim,8,"");
        rhsState.AddField("scalarGrad",'n',numDim*numScalar,8,"");
        tempVector = new double [numPointsBuffer*numDim];
        tempScalar = new double [numPointsBuffer];
      }
      rhsState.Create(numPointsBuffer,0);
      velHat = rhsState.template GetFieldBuffer<double>("vHat");
      viscidEnergy = rhsState.template GetFieldBuffer<double>("viscidEnergy");
      if(artificialDissipation){
        if(velDiv){
          rhsState.SetFieldBuffer("velDiv",velDiv);
        } else {
          velDiv           = rhsState.template GetFieldBuffer<double>("velDiv");
        }
        sigmaDissipation = rhsState.template GetFieldBuffer<double>("sigmaDiss");
        alphaDissipation = rhsState.template GetFieldBuffer<double>("alphaDiss");
      } else {
        sigmaDissipation = NULL;
        alphaDissipation = NULL;
      }

      velocityHandle = inState.GetDataIndex("velocity");
      if(velocityHandle >= 0){
        pcpp::field::dataset::DataBufferType &velocityData(inState.Field(velocityHandle));
        velocityPtr = velocityData.Data<double>();
      }

      // inviscid flux buffers
      if(inviscidFluxBuffers.empty()){
        inviscidFluxBuffers.resize(numEquations);
      }
      int inviscidFluxHandle = rhsState.GetDataIndex("inviscidFlux");
      inviscidFluxBuffer = rhsState.GetRealFieldBuffer(inviscidFluxHandle);
      for(int i=0;i<numEquations;i++){
        inviscidFluxBuffers[i] = &inviscidFluxBuffer[i*numPointsBuffer];
      }

      // viscid flux buffers
      if(viscidFluxBuffers.empty()){
        viscidFluxBuffers.resize(numEquations);
      }
      int viscidFluxHandle = rhsState.GetDataIndex("viscidFlux");
      viscidFluxBuffer = rhsState.GetRealFieldBuffer(viscidFluxHandle);
      for(int i=0;i<numEquations;i++){
        viscidFluxBuffers[i] = &viscidFluxBuffer[i*numPointsBuffer];
      }

      // flux deriviatives
      if(dFluxBuffers.empty()){
        dFluxBuffers.resize(numEquations);
      }
      int dFluxHandle = rhsState.GetDataIndex("dFlux");
      dFluxBuffer = rhsState.GetRealFieldBuffer(dFluxHandle);
      for(int i=0;i<numEquations;i++){
        dFluxBuffers[i] = &dFluxBuffer[i*numPointsBuffer];
      }

      if(refRe > 0){
        if(velGradBuffers.empty()){
          velGradBuffers.resize(numDim);
        }
        int velGradHandle = rhsState.GetDataIndex("velGrad");
        velGradBufferPtr = rhsState.GetRealFieldBuffer(velGradHandle);
        for(int i=0;i<numDim;i++){
          velGradBuffers[i] = &velGradBufferPtr[i*numDim*numPointsBuffer];
        }
        if(numScalar > 0){
          if(scalarGrad == NULL)
            scalarGrad = rhsState.GetRealFieldBuffer(rhsState.GetDataIndex("scalarGrad"));
          else
            rhsState.SetFieldBuffer("scalarGrad",scalarGrad);
        }
        if(tauBuffers.empty()){
          tauBuffers.resize(numDim*numDim);
        }
        int tauHandle = rhsState.GetDataIndex("tau");
        tauBufferPtr = rhsState.GetRealFieldBuffer(tauHandle);
        
        // tau is a symmetric 2-tensor, only store the upper half
        // and set the pointers of the lower half to the memory 
        // for the upper half
        size_t bufferIndex = 0;
        for(int i=0;i<numDim;i++){
          for(int j=0;j<numDim;j++){
            int index = i*numDim + j; 
            if (i <= j){
              tauBuffers[index] = &tauBufferPtr[bufferIndex];
              bufferIndex += numPointsBuffer;
            } else {
              tauBuffers[index] = tauBuffers[j*numDim+i];
            }
          }
        }

        // thermal conduction storage
        if(temperatureGradBuffers.empty()){
          temperatureGradBuffers.resize(numDim);
        }
        int temperatureGradHandle = rhsState.GetDataIndex("temperatureGrad");
        temperatureGradBufferPtr = rhsState.GetRealFieldBuffer(temperatureGradHandle);
        for(int i=0;i<numDim;i++){
          temperatureGradBuffers[i] = &temperatureGradBufferPtr[i*numPointsBuffer];
        }
 
        if(heatFluxBuffers.empty()){
          heatFluxBuffers.resize(numDim);
        }
        int heatFluxHandle = rhsState.GetDataIndex("heatFlux");
        heatFluxBufferPtr = rhsState.GetRealFieldBuffer(heatFluxHandle);
        
        for(int i=0;i<numDim;i++){
          heatFluxBuffers[i] = &heatFluxBufferPtr[i*numPointsBuffer];
        }
  
      }

      // 
      // register needed buffers with the eos object
      //
      eosPtr->SetupPressureBuffer(pressurePtr);
      eosPtr->SetupTemperatureBuffer(temperaturePtr);
      eosPtr->SetupSpecificVolumeBuffer(rhom1Ptr);
      eosPtr->SetupInternalEnergyBuffer(internalEnergyPtr);

      dScalarHandle = scalarHandle;

      stateInitialized = true;

      return(0);
    };
    
    // *NOT* thread-safe, call on master only
    int SetTargetState(StateType &inState)
    {
      if(!stateInitialized)
        return(1);

      targetState.CopyStateData(inState);

      int rhoTargetHandle = inState.GetDataIndex("rhoTarget");
      int rhoETargetHandle   = -1;
      int rhoVTargetHandle   = -1;
      int scalarTargetHandle = -1;

      bool targetData = false;
      // If no acutal target data, then use the input state as the target
      if(rhoTargetHandle < 0){
        rhoTargetHandle  = targetState.GetDataIndex("rho");
        rhoVTargetHandle = targetState.GetDataIndex("rhoV");
        rhoETargetHandle = targetState.GetDataIndex("rhoE");
        if(numScalar > 0){
          scalarTargetHandle = targetState.GetDataIndex("scalarVars");
        }
      } else {
        targetData = true;
        rhoVTargetHandle = inState.GetDataIndex("rhoVTarget");
        rhoETargetHandle = inState.GetDataIndex("rhoETarget");
        if(numScalar > 0){
          scalarTargetHandle = inState.GetDataIndex("scalarVarsTarget");
        }
      }
      if(rhoTargetHandle  < 0   ||
         rhoVTargetHandle < 0   ||
         rhoETargetHandle < 0   || 
         (numScalar > 0 && scalarTargetHandle < 0)){
        std::cerr << "RHS Target Handles problem: " << rhoTargetHandle << " "
                  << rhoVTargetHandle << " " << rhoETargetHandle << std::endl;
        return(1);
      }
      
      // Set up state field buffers if necessary
      if(targetData){
        targetState.SetFieldBuffer("rho",inState.GetRealFieldBuffer(rhoTargetHandle));
        targetState.SetFieldBuffer("rhoV",inState.GetRealFieldBuffer(rhoVTargetHandle));
        targetState.SetFieldBuffer("rhoE",inState.GetRealFieldBuffer(rhoETargetHandle));
        if(numScalar > 0)
          targetState.SetFieldBuffer("scalarVars",inState.GetRealFieldBuffer(scalarTargetHandle));
      }
      
      rhoTargetPtr         = targetState.GetRealFieldBufferByHandle(rhoHandle);
      rhoVTargetPtr        = targetState.GetRealFieldBufferByHandle(rhoVHandle);
      rhoETargetPtr        = targetState.GetRealFieldBufferByHandle(rhoEHandle);
      if(numScalar > 0){
        scalarTargetPtr = targetState.GetRealFieldBufferByHandle(scalarHandle);
      }

      targetStateBuffers.resize(0);
      int numTargetVars = numEquations+1;
      targetStateBuffers.resize(numTargetVars,NULL);

      int iEquation = 0;
      targetStateBuffers[iEquation++] = rhoTargetPtr;
      for(int iDim = 0;iDim < numDim;iDim++)
        targetStateBuffers[iEquation++] = &rhoVTargetPtr[iDim*numPointsBuffer];
      targetStateBuffers[iEquation++] = rhoETargetPtr;
    
      if(numScalar > 0){
        for(int iScalar = 0;iScalar < numScalar;iScalar++)
          targetStateBuffers[iEquation++] = &scalarTargetPtr[iScalar*numPointsBuffer];
      }

      return(0);

    };
    
    
    int SetState(StateType &inState)
    {
      if(!stateInitialized)
        return(1);
      
      rhoPtr         = inState.GetRealFieldBufferByHandle(rhoHandle);
      rhoVPtr        = inState.GetRealFieldBufferByHandle(rhoVHandle);
      rhoEPtr        = inState.GetRealFieldBufferByHandle(rhoEHandle);

      int iEquation = 0;
      myStateBuffers[iEquation++] = rhoPtr;
      for(int iDim = 0;iDim < numDim;iDim++)
        myStateBuffers[iEquation++] = &rhoVPtr[iDim*numPointsBuffer];
      myStateBuffers[iEquation++] = rhoEPtr;
      
      if(numScalar > 0){
        scalarPtr = inState.GetRealFieldBufferByHandle(scalarHandle);
        for(int iScalar = 0;iScalar < numScalar;iScalar++)
          myStateBuffers[iEquation++] = &scalarPtr[iScalar*numPointsBuffer];
      }
      
      
      //       timePtr        = inState.GetRealFieldBuffer(timeHandle);
      //       time = *timePtr;
      
      //       if(dtHandle >= 0){
      //         dtPtr  = inState.GetRealFieldBuffer(dtHandle);
      //       }
      //       if(cflHandle >= 0){
      //         cflPtr = inState.GetRealFieldBuffer(cflHandle);
      //       }
      
      return(0);

    };
    
    int SetRHSState(StateType &rhsState)
    {

      if(dRhoHandle < 0){
        dRhoHandle  = rhsState.GetFieldHandle("rho");
        dRhoVHandle = rhsState.GetFieldHandle("rhoV");
        dRhoEHandle = rhsState.GetFieldHandle("rhoE");
        if(numScalar > 0)
          dScalarHandle = rhsState.GetFieldHandle("scalarVars");
      }
      
      rhoRHSPtr         = rhsState.GetRealFieldBufferByHandle(dRhoHandle);
      rhoVRHSPtr        = rhsState.GetRealFieldBufferByHandle(dRhoVHandle);
      rhoERHSPtr        = rhsState.GetRealFieldBufferByHandle(dRhoEHandle);

      if(numScalar > 0)
        scalarRHSPtr = rhsState.GetRealFieldBufferByHandle(dScalarHandle);
      else
        scalarRHSPtr = NULL;

      int iEquation = 0;
      myRHSBuffers[iEquation++] = rhoRHSPtr;
      for(int iDim = 0;iDim < numDim;iDim++)
        myRHSBuffers[iEquation++] = &rhoVRHSPtr[numPointsBuffer*iDim];
      myRHSBuffers[iEquation++] = rhoERHSPtr;
      for(int iScalar = 0;iScalar < numScalar;iScalar++)
        myRHSBuffers[iEquation++] = &scalarRHSPtr[numPointsBuffer*iScalar];
      
      return(0);    

    };

    //  ============= Stencil & Operator utils ========

    // This routine *needs* to be thread-aware - but calling sequence needs to be 
    // worked out.  Currently is called before threads are initialized.
    int SetupStencilConnectivities()
    {

      if(!operatorInitialized)
        return(1);

      int myThreadId = 0;
      int numThreads = 1;

      // #ifdef USE_OMP
      //       myThreadId = omp_get_thread_num();
      //       numThreads = omp_get_num_threads();
      //       masterThread = (myThreadId == 0);
      // #endif
      
      //      if(masterThread)

#pragma omp master
      {
        if(!stencilConn){
          stencilConn = new int [numDim*numPointsBuffer];
        }
        if(artificialDissipation){
          artDissConn      = new int [numDim*numPointsBuffer];
          artDissConnSplit = new int [numDim*numPointsBuffer];
        }
        if(!iMask){
          iMask = new int [numPointsBuffer];
        }
      }
          
      //#pragma omp barrier
      
      // Create stencil connectivity
      //      if(numThreads == 1){
      // Only deals with actual domain boundaries - no holes
      plascom2::operators::sbp::CreateStencilConnectivity(numDim,&bufferSizes[0],&opInterval[0],
                                                          boundaryDepth,&periodicDirs[0],
                                                          stencilConn,true);

      if(artificialDissipation){
        plascom2::operators::sbp::CreateStencilConnectivity(numDim,&bufferSizes[0],&opInterval[0],
                                                            artDissOperator.boundaryDepth,&periodicDirs[0],
                                                            artDissConn,true);

        plascom2::operators::sbp::CreateStencilConnectivity(numDim,&bufferSizes[0],&opInterval[0],
                                                            artDissOperatorSplit.boundaryDepth,&periodicDirs[0],
                                                            artDissConnSplit,true);
      }
      //      } else {
      // plascom2::operators::sbp::CreateStencilConnectivity(numDim,&bufferSizes[0],
      //                         &fortranPartitionBufferIntervals[myThreadId][0],
      //                         boundaryDepth,&periodicDirs[0],
      //                         stencilConn,true); 
      //    }
    
    
    // plascom2::operators::sbp::CreateStencilConnectivity(numDim,&gridSize[0],opInterval,boundaryDepth,stencilConn,true);
    
    // bug! unlike stencilConn, these need to be thread-specific
    // @todo Thread-specific dual stencil connectivities in RHS
    // Invert the connectivity
    // plascom2::operators::InvertStencilConnectivity(numDim,opInterval,stencilConn,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)]...]
    
    //       plascom2::operators::sbp::InvertStencilConnectivity(numDim,&bufferSizes[0],&opInterval[0],numStencils,
    //                                                           stencilConn,dualStencilConn,numPointsStencil,true);
      return(0);
      
    };
    

    int SetupOperators(OperatorType &inOperator)
    {

      int myThreadId = 0;
      int numThreads = 1;
      
      // #ifdef USE_OMP
      //       myThreadId = omp_get_thread_num();
      //       masterThread = (myThreadId == 0);
      // #endif
      int errorCode = 0;

#pragma omp master
      {
        // ==== Main differential operators ===
        if(!operatorInitialized){
          myOperator.Copy(inOperator);
          
          numStencils      = myOperator.numStencils;
          stencilStarts    = myOperator.stencilStarts; 
          stencilOffsets   = myOperator.stencilOffsets;
          stencilSizes     = myOperator.stencilSizes;
          stencilWeights   = myOperator.stencilWeights;
          boundaryDepth    = myOperator.boundaryDepth;
          boundaryWeight   = myOperator.boundaryWeight;
          numStencilValues = myOperator.numValues;
          
          // Forticate the stencil starts
          for(int iStencil = 0;iStencil < numStencils;iStencil++)
            stencilStarts[iStencil]++;
          operatorInitialized = true;
          // === artificial dissipation operators ===
          if(artificialDissipation){
            int interiorOrderAD = (myOperator.overallOrder-1)*2;
            if(plascom2::operators::dissipation::Initialize(artDissOperator,artDissOperatorSplit,interiorOrderAD)){
              std::cerr << "RHS ERROR: Artificial dissipation operator failed to initialize." << std::endl;
              errorCode = 1;
              operatorInitialized = false;
            } else {
              int *startPtr = artDissOperator.stencilStarts;
              int numStencilsAD = artDissOperator.numStencils;
              // Forticate the stencil starts
              for(int iStencil = 0;iStencil < numStencilsAD;iStencil++)
                startPtr[iStencil]++;
              startPtr = artDissOperatorSplit.stencilStarts;
              numStencilsAD = artDissOperator.numStencils;
              for(int iStencil = 0;iStencil < numStencilsAD;iStencil++)
                startPtr[iStencil]++;
            }
          } 
        }
      }

      if(errorCode > 0){
        operatorInitialized = false;
        return(1);
      }

      //#pragma omp barrier
      if(SetupStencilConnectivities())
        return(1);
      
      gridPtr->SetupDifferentialOperator(&myOperator,&stencilConn[0]);
      std::ostringstream metricMessages;
      if(gridPtr->ComputeMetrics(metricMessages)){
        std::ostringstream failMessage;
        failMessage << "euler::rhs::SetupOperators: Error: grid::ComputeMetrics failed "
                    << "with messages:" << std::endl << metricMessages.str() << std::endl;
        if(globalPtr){
          globalPtr->ErrOut(failMessage.str());
        } else {
          std::cout << failMessage.str();
        }
        return(1);
      } else {
        if(globalPtr){
          globalPtr->StdOut("Grid metric messages:\n",3);
          globalPtr->StdOut(metricMessages.str());
        } else {
          std::cout << "Grid metric messages:\n" << metricMessages.str() << std::endl;
        }
      }
      if(baseStatePtr){
        int gridMetricsHandle = baseStatePtr->GetDataIndex("gridMetrics");
        if(gridMetricsHandle >= 0){
          std::vector<double> gridMetric(gridPtr->Metric());
          baseStatePtr->SetFieldBuffer("gridMetrics",&gridMetric[0]);
        }
      }
            
      return(0);
      
    };
    
    // ====== Thread set up/handling =====

    int InitThreadIntervals()
    {
      int myThreadId = 0;
      int numThreads = 1;
      bool inParallel = false;

#ifdef USE_OMP
      inParallel = omp_in_parallel();
      if(inParallel){
        myThreadId = omp_get_thread_num();
        numThreads = omp_get_num_threads();
      }
#endif
      
#pragma omp master
      {
        fortranBufferIntervals.resize(numThreads);
        fortranPartitionBufferIntervals.resize(numThreads);
      }
      
      if(inParallel){
#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);

    };
    

    //  ====== Grid (and halo) handling =====
    int SetGrid(GridType &inGrid)
    {
      
      int numThreads = 1;

#ifdef USE_OMP
      numThreads = omp_get_num_threads();
#endif

      if(gridPtr != &inGrid){

        int myThreadId = 0;

#ifdef USE_OMP
        myThreadId  = omp_get_thread_num();
#endif
      

#pragma omp master
        {

          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);
          gridType                = gridPtr->Type();
          
          //           std::vector<double> &gridJacobianVector(gridPtr->Jacobian());
          //           std::vector<double> &gridMetricVector(gridPtr->Metric());
          
          //           gridJacobian = &gridJacobianVector[0];
          //           gridMetric   = &gridMetricVector[0];
          //           if(velHat)
          //             delete [] velHat;
          //           velHat = new double [numDim*numPointsBuffer];
          gridInitialized    = true;
          //           gridDX.resize(numDim,1.0);
          //           gridDXM1.resize(numDim,1.0);
          //           gridMetric.resize(numDim,1.0);
          //           gridJacobian = 1.0;
          //           std::vector<double> &dX(inGrid.DX());
          //           std::vector<double>::iterator dXIt = dX.begin();
          //           std::vector<double>::iterator gridDXIt = gridDX.begin();
          //           std::vector<double>::iterator gridDXM1It = gridDXM1.begin();
          //           while(dXIt != dX.end()){
          //             double dx = *dXIt++;
          //             *gridDXIt++   = dx;
          //             *gridDXM1It++ = 1/dx;
          //             gridJacobian *= 1/dx;
          //           }
          
          //           dXIt = dX.begin();
          //           gridDXM1It = gridMetric.begin();
          //           while(dXIt != dX.end()){
          //             double dx = *dXIt++;
          //             *gridDXM1It++ = 1.0/(dx*gridJacobian);
          //           }
          
          partitionBufferInterval.Flatten(opInterval);
          
          // 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();
      }
    
      //       else if ((numThreads > 1) && !threadsInitialized) { // JIC not already done
      //         InitThreadIntervals();
      //       }
      return(0);
    };


    int HaloSetup(HaloType &inHalo)
    {
      
#pragma omp master
      {
        SetHalo(inHalo);
        periodicDirs = inHalo.PeriodicDirs();
        if(periodicDirs.empty()){
          periodicDirs.resize(numDim,1);
        }
        //         int pDim = periodicDirs.size();
        //         std::cout << "Periodic(";
        //         for(int iDim = 0;iDim < pDim;iDim++){
        //           if(iDim > 0) std::cout << ",";
        //           std::cout << periodicDirs[iDim];
        //         }
        //         std::cout << ")" << std::endl;

        maskMessageId = inHalo.CreateMessageBuffers(1,4);
        if(exchangeDV)
          dvMessageId   = inHalo.CreateMessageBuffers(numDV);
        if(exchangeFlux)
          fluxMessageId = inHalo.CreateMessageBuffers(numEquations);
        scalarMessageId = inHalo.CreateMessageBuffers(1);
      }
      return(0);      
    };
  
    int SetHalo(HaloType &inHalo)
    {
      haloPtr = &inHalo;
      return(0);
    };


    // ===== Boundary related stuff ====

    // 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.
    // Threaded fail 
    void SetDomainBCs(std::vector<BoundaryType> &domainBoundaries,std::vector<BCType> &domainBCs)
    { 
      
      domainBCsPtr = &domainBCs;
      domainBoundariesPtr = &domainBoundaries;      
      assert(subRegionsPtr != NULL);
    
      std::vector<GridRegionType> &gridRegions(*subRegionsPtr);

      int numBCs               = domainBCs.size();
      int numGridRegions       = gridRegions.size();
      int numDomainBoundaries  = domainBoundaries.size();
      int holeStencil          = myOperator.numStencils;
      int adHole               = artDissOperator.numStencils;
      int adSplitHole          = artDissOperatorSplit.numStencils;
      int adBoundaryDepth      = artDissOperator.boundaryDepth;
      int adSplitBoundaryDepth = artDissOperatorSplit.boundaryDepth;

      //       std::cout << "Boundary stencil processing report: "  << std::endl
      //                 << "HoleStencil:      " << holeStencil     << std::endl
      //                 << "ADHoleStencil:    " << adHole          << std::endl
      //                 << "ADSplitHole:      " << adSplitHole     << std::endl
      //                 << "ADBoundaryDepth:  " << adBoundaryDepth << std::endl
      //                 << "ADBoundaryDepthS: " << adSplitBoundaryDepth << std::endl;
      
      //       if(numBCs == 0)
      //         return;
      //       if(numDomainBoundaries == 0)
      //         return;
      
      
      pcpp::IndexIntervalType bufferInterval(numDim,&bufferSizes[0]);
      std::vector<size_t> partitionBufferExtent;
      partitionBufferInterval.Flatten(partitionBufferExtent);
      
      
      // Initialize the mask to zero
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
        iMask[iPoint] = 0;
      bool haveHoles = false;
      
      // Check for legacy IBLANK (0 = hole)
      size_t numHoles = gridPtr->FlagMask(iMask,0,holeMask);
      size_t numHolesGrid = 0;
      gridCommPtr->AllReduce(numHolesGrid,numHoles,ix::comm::DTSIZET,ix::comm::SUMOP);
      if(numHolesGrid > 0){
        haveHoles = true;
      }

      // First, loop through all holes and create the iMask
      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;
        BCType &boundaryCondition(domainBoundary.BC());
        // unique int id for bc
        int bcType = boundaryCondition.BCType();
        
        if(bcType == navierstokes::HOLE) { 
          // set hole bits in mask
          haveHoles = true;
          pcpp::IndexIntervalType opInterval;
          boundaryPartitionInterval.RelativeTranslation(partitionInterval,
                                                        partitionBufferInterval,
                                                        opInterval);
          if(opInterval.NNodes() > 0) {
            std::vector<size_t> opExtents;
            opInterval.Flatten(opExtents);
            mask::SetMask(bufferSizes,opInterval,holeMask,iMask);
          }
        }
      }
      
      // MJA
      // Bug here? Only needed for artificial dissipation
      //
      //      std::cout << "passed this part" << std::endl;
      if(haveHoles){

        int returnCode = 0;

        haloPtr->PackMessageBuffers(maskMessageId,0,1,iMask);
        haloPtr->SendMessage(maskMessageId,gridNeighbors,*gridCommPtr);
        haloPtr->ReceiveMessage(maskMessageId,gridNeighbors,*gridCommPtr);
        haloPtr->UnPackMessageBuffers(maskMessageId,0,1,iMask);
        
        returnCode = plascom2::operators::sbp::DetectHoles(numDim,&bufferSizes[0],&partitionBufferExtent[0],
                                                           boundaryDepth,holeMask,iMask,stencilConn);
        
            
        if(artificialDissipation){
          returnCode = plascom2::operators::sbp::DetectHoles(numDim,&bufferSizes[0],&partitionBufferExtent[0],
                                                             adBoundaryDepth,holeMask,iMask,artDissConn);
          returnCode = plascom2::operators::sbp::DetectHoles(numDim,&bufferSizes[0],&partitionBufferExtent[0],
                                                             adSplitBoundaryDepth,holeMask,iMask,artDissConnSplit);
        }       
        if(returnCode){
          std::string errorMessage("rhs::SetDomainBCs:Error: Hole detection failed.\n");
          if(globalPtr){
            globalPtr->ErrOut(errorMessage);
          } else {
            std::cerr << errorMessage;
          }
        }
      }

    };
      
      //  == auxiliary methods ==
    int InitStep(GridType &inGrid,StateType &inState)
    {
      return(0);
    };    
    virtual int ProcessInterior2(GridType &inGrid,StateType &inState,StateType &rhsState)
    { return(0); };
    virtual int ProcessHalo(HaloType &inHalo,GridType &inGrid,StateType &inState,StateType &rhsState)
    { return(0); };
    virtual int ProcessInterior(GridType &inGrid,StateType &inState,StateType &rhsState)
    { return(0); };
    virtual int ProcessBoundary(GridType &inGrid,StateType &inState,StateType &rhsState)
    { return(0); };

    std::vector<double *> &DVBuffers()
    { return(dvBuffers); };
    std::vector<double *> &TVBuffers()
    { return(tvBuffers); };
    std::vector<double *> &StateBuffers()
    { return(myStateBuffers); };
    std::vector<double *> &RHSBuffers()
    { return(myRHSBuffers); };
    double *InviscidFluxBuffer()
    { return(inviscidFluxBuffer); };
    double *DFluxBuffer()
    { return(dFluxBuffer); };
    int NumberOfEquations()
    { return(numEquations); };

    
    WENO::CoeffsWENO &WENOCoefficients()
    { return(coeffsWENO); };
    std::vector<double *> &VelGradBuffers()
    { return(velGradBuffers); };
    std::vector<double *> &TemperatureGradBuffers()
    { return(temperatureGradBuffers); };
    
    void InitMinSpacing()
    {
#pragma omp master
      {
        minSpacing.resize(numDim,std::numeric_limits<double>::max());
        //        if(gridType < simulation::grid::RECTILINEAR){
        if(true){
          for(int iDim = 0;iDim < numDim;iDim++){
            minSpacing[iDim] = gridDX[iDim];
          }
          
          //         } else {
          //           std::vector<double> &coordinateData(gridPtr->CoordinateData());
          //           std::vector<size_t> dimOffset(numDim,0);
          //           for(int iDim = 0;iDim < numDim;iDim++)
          //             dimOffset[iDim] = iDim * numPointsBuffer;
          
          //           pcpp::IndexIntervalType globalInterval;
          //           globalInterval.InitSimple(gridSizes);
          //           std::vector<size_t> partitionPointIndices;
          //           globalInterval.GetFlatIndices(partitionBufferInterval,partitionPointIndices);
          //           std::vector<size_t>::iterator pointIt = partitionPointIndices.begin();
          //           while(pointIt != partitionPointIndices.end()){
          //             size_t bufferIndex = *pointIt++;
          //             for(int iDim = 0;iDim < numDim;iDim++){
          //               if(coordinateData[bufferIndex+dimOffset[iDim]] < minSpacing[iDim])
          //                 minSpacing[iDim] = coordinateData[bufferIndex+dimOffset[iDim]];
          //             }
          //           }
        }    
      }
    };
    
    // Must be called on all threads
    double TimeStep(double *outCFL = NULL)
    {
      //       if(minSpacing.empty())
      //         InitMinSpacing();
      //       double spaceFactor = 0.0;
      //       for(int iDim = 0;iDim < numDim;iDim++)
      //         spaceFactor += 1.0/(minSpacing[iDim]);
      
      *dtPtr = 0;
      if(timeSteppingMode == CONSTANT_DT){
#pragma omp master
        {
          if(*dtPtr != inputDT){
            *dtPtr  = inputDT;
            *cflPtr = *dtPtr;
            //            *cflPtr = refSndSpd*(*dtPtr)*spaceFactor;
          }
        }
#pragma omp barrier
        if(outCFL != NULL){
          *outCFL = *cflPtr;
        }
      } 
      //       else if (timeSteppingMode == CONSTANT_CFL){
      
      // #pragma omp master
      //         {
      //           if(*cflPtr != inputCFL){
      //             *cflPtr = inputCFL;
      //             *dtPtr = *cflPtr/(refSndSpd*spaceFactor);
      //           }
      //         }
      // #pragma omp barrier
      
      //       }
      return(*dtPtr);
    };
    
  private:
    int Init()
    {

      int myThreadId = 0;
      int numThreads = 1;

#ifdef USE_OMP
      myThreadId = omp_get_thread_num();
      numThreads = omp_get_num_threads();
#endif

#pragma omp master 
      {
        myCFL              = 1.0;
        myDT               = 0.0;

        //      if(masterThread){
        // Conserved state 
        rhoHandle          = -1;
        rhoVHandle         = -1;
        rhoEHandle         = -1;
        scalarHandle       = -1;
        dRhoHandle         = -1;
        dRhoVHandle        = -1;
        dRhoEHandle        = -1;
        dScalarHandle      = -1;
        
        // Dependent state
        temperatureHandle  = -1;
        pressureHandle     = -1;
        gradVHandle        = -1;
        gradTHandle        = -1;
        
        // Parameters
        gammaHandle        = -1;
        cflHandle          = -1;
        dtHandle           = -1;
        numScalar          =  0;
        sigmaHandle        = -1;
        cHandle            = -1;
        reHandle           = -1;
        cpHandle           = -1;
        gasModel           =  0;
        transportModel     =  0;
        gridType           =  0;

        // Optional data handles
        velocityHandle     = -1;
        rhom1Handle        = -1;
        internalEnergyHandle = -1;

        // Misc
        numPointsPartition = 0;
        numPointsBuffer    = 0;
        numDim             = 0;
        numStateFields     = 3;
        numDV              = 6;
        numStateVar        = 5;
        baseStatePtr       = NULL;
        gridPtr            = NULL;
        cPtr               = NULL;
        rePtr              = NULL;
        cpPtr              = NULL;
        timePtr            = NULL;
        rhoPtr             = NULL;
        rhoVPtr            = NULL;
        rhoEPtr            = NULL;
        scalarPtr          = NULL;
        rhoTargetPtr       = NULL;
        rhoVTargetPtr      = NULL;
        rhoETargetPtr      = NULL;
        scalarTargetPtr    = NULL;
        pressurePtr        = NULL;
        temperaturePtr     = NULL;
        rhoRHSPtr          = NULL;
        rhoVRHSPtr         = NULL;
        rhoERHSPtr         = NULL;
        scalarRHSPtr       = NULL;
        massFluxPtr        = NULL;
        momentumFluxPtr    = NULL;
        energyFluxPtr      = NULL;
        scalarFluxPtr      = NULL;
        fluxComponentPtr   = NULL;
        velocityPtr        = NULL;
        rhom1Ptr           = NULL;
        internalEnergyPtr  = NULL;
        massSourcePtr      = NULL;
        momentumSourcePtr  = NULL;
        energySourcePtr    = NULL;
        scalarSourcePtr    = NULL;
        dvBuffer           = NULL;
        viscidFluxBuffer   = NULL;
        inviscidFluxBuffer = NULL;
        dFluxBuffer        = NULL;
        //        scaledPressure     = NULL;
        scaledTau          = NULL;
        velHat             = NULL;
        viscidEnergy       = NULL;
        velGradBufferPtr   = NULL;
        velDiv             = NULL;
        sigmaDissipation   = NULL;
        alphaDissipation   = NULL;
        scalarGrad         = NULL;
        tempVector         = NULL;
        tempScalar         = NULL;
        inputDiffusivity   = NULL;
        scalarDiffusivity  = NULL;
        // EOS
        eosPtr             = NULL;

        temperatureGradBufferPtr   = NULL;
        cflPtr             = &myCFL;
        dtPtr              = &myDT;
        timeSteppingMode   = CONSTANT_DT;

        // BCs-related
        subRegionsPtr       = NULL;
        domainBoundariesPtr = NULL;
        domainBCsPtr        = NULL;

        // Operator-related
        stencilConn         = NULL;
        artDissConn         = NULL;
        artDissConnSplit    = NULL;
        iMask               = NULL;

        //     dvBuffers           = NULL;
        //     fluxBuffers         = NULL;
        stateInitialized    = false;
        gridInitialized     = false;
        threadsInitialized  = false;
        operatorInitialized = false;
        errorCode = 0;
        
        dvMessageId     = -1;
        fluxMessageId   = -1;
        maskMessageId   = -1;
        scalarMessageId = -1;

      }
      return 0;
    };

    //   public:
    //     ~rhs()
    //     {
    //       std::cout << "Enter destructor" << std::endl;
    //       if(fluxBuffer)
    //         delete [] fluxBuffer;
    //       if(dFluxBuffer)
    //         delete [] dFluxBuffer;
    //       if(scaledPressure)
    //         delete [] scaledPressure;
    //       for(int iDV = 0;iDV < numDV;iDV++)
    //         if(dvBuffersMine[iDV])
    //           delete [] dvBuffers[iDV];
    //       if(velHat)
    //         delete [] velHat;
    //       std::cout << "Exit destructor" << std::endl;
    //     };
    
  protected:

    // State ptr
    StateType *baseStatePtr;

    // Data handles
    int rhoHandle;
    int rhoVHandle;
    int rhoEHandle;
    int dRhoHandle;
    int dRhoVHandle;
    int dRhoEHandle;
    int scalarHandle;
    int dScalarHandle;
    int timeHandle;
    int gammaHandle;
    int powerHandle;
    int betaHandle;
    int bulkViscFacHandle;
    int sigmaHandle;
    int cflHandle;
    int dtHandle;
    int pressureHandle;
    int temperatureHandle;
    int muHandle;
    int lambdaHandle;
    int kappaHandle;
    int rhom1Handle;
    int internalEnergyHandle;
    int velocityHandle;
    int gradVHandle;
    int gradTHandle;
    int cHandle;
    int reHandle;
    int cpHandle;

    // Parameters
    double time;
    double myDT;
    double myCFL;
    double inputCFL;
    double inputDT;

    int nonDimensionalMode;
    // transport model specific
    double power;
    double beta;
    double bulkViscFac;
    double refTemperature;
    double refPressure;
    double refRho;
    double refEnergy;
    double cfl;
    double refSndSpd;
    double refVelocity;
    double refLength;
    double refCp;
    double refCv;
    double refRe;
    double refPr;
    double refMu;
    double refKappa;
    double refSpecGasConst;
    double pressureScale;
    double temperatureScale;
    double velocityScale;
    double specGasConst;
    double inputSigma;
    double inputSigmaDil;

    int    gasModel;
    int    transportModel;
    const double *powerPtr;
    const double *betaPtr;
    const double *bulkViscFacPtr;
    const double *timePtr;
    //    const double *cflPtr;
    const double *sigmaPtr;
    const double *sigmaDilPtr;
    const double *cPtr;
    const double *rePtr;
    const double *prPtr;
    const double *cpPtr;
    const double *inputDiffusivity;

    int numScalar;
    int numStateFields;
    int numStateVar;
    int numDV;
    int numTV;
    int timeSteppingMode;

    // Grid data
    GridType *gridPtr;
    //     StateType *targetStatePtr;
    //     StateType targetState;
    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> gridDX;
    std::vector<double> gridDXM1;
    std::vector<double> minSpacing;
    int gridType;

    double gridJacobian;
    std::vector<double> gridMetric;

    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;
    StateType targetState;
    StateType rhsState;


    // input State information
    double *rhoPtr;
    double *rhoVPtr;
    double *rhoEPtr;
    double *scalarPtr;

    double *rhoTargetPtr;
    double *rhoVTargetPtr;
    double *rhoETargetPtr;
    double *scalarTargetPtr;
    double *pressurePtr;
    double *temperaturePtr;
    double *muPtr;
    double *lambdaPtr;
    double *kappaPtr;
    double *scalarDiffusivity;

    double *dtPtr;
    double *cflPtr;

    // rhs State information
    double *rhoRHSPtr;
    double *rhoVRHSPtr;
    double *rhoERHSPtr;
    double *scalarRHSPtr;
  
    // flux storage
    double *massFluxPtr;
    double *momentumFluxPtr;
    double *energyFluxPtr;
    double *scalarFluxPtr;
    double *fluxComponentPtr;

    // Auxiliary data
    double *velocityPtr;
    double *rhom1Ptr;
    double *internalEnergyPtr;
    double *velHat;
    double *viscidEnergy;
    double *velDiv;
    double *sigmaDissipation;
    double *alphaDissipation;

    // Source data
    double *massSourcePtr;
    double *momentumSourcePtr;
    double *energySourcePtr;
    double *scalarSourcePtr;
  
    // Halo data storage
    std::vector< std::vector<double *> > haloDVBuffers;
    std::vector< std::vector<double *> > haloFluxBuffers;
    std::vector<std::vector<double *> > haloStateBuffers;

    // My storage
    std::vector<bool> dvBuffersMine;
    std::vector<double *> dvBuffers;
    std::vector<double *> tvBuffers;
    std::vector<double *> inviscidFluxBuffers;
    std::vector<double *> viscidFluxBuffers;
    std::vector<double *> dFluxBuffers;
    std::vector<double *> myStateBuffers;
    std::vector<double *> targetStateBuffers;
    std::vector<double *> myRHSBuffers;
    std::vector<double *> velGradBuffers; // dv_i/dx_j
    std::vector<double *> temperatureGradBuffers; // dT/dx_i
    std::vector<double *> tauBuffers; 
    std::vector<double *> heatFluxBuffers; 
    double *dvBuffer;
    double *tvBuffer;
    double *inviscidFluxBuffer;
    double *viscidFluxBuffer;
    double *dFluxBuffer;
    //    double *scaledPressure;
    double *scaledTau;
    double *velGradBufferPtr;
    double *temperatureGradBufferPtr;
    double *tauBufferPtr;
    double *heatFluxBufferPtr;
    double *scalarGrad;
    double *tempVector;
    double *tempScalar;

    // Operator
    OperatorType myOperator;
    OperatorType artDissOperator;
    OperatorType artDissOperatorSplit;

    int numStencils;
    int boundaryDepth;
    int numStencilValues;
    int *stencilSizes;
    int *stencilOffsets;
    int *stencilConn;
    int *artDissConn;
    int *artDissConnSplit;
    int *iMask;
    int *stencilStarts;
    double *stencilWeights;
    size_t *dualStencilConn;
    size_t *numPointsStencil;

    // bookkeeping
    bool stateInitialized;
    bool gridInitialized;
    bool threadsInitialized;
    bool operatorInitialized;
    std::vector<bool> haveHalo;
    int errorCode;
    
    // message Ids
    int dvMessageId;
    int fluxMessageId;
    int maskMessageId;
    int scalarMessageId;
    //    int numDVMessageComponents;
    int numEquations;

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

    // Processing options
    //    bool useAlternateRHS;
    bool exchangeFlux;
    bool artificialDissipation;
    bool exchangeDV;
    bool useWENO;
   
    // WENO data structures
    WENO::CoeffsWENO coeffsWENO;

    // equation of state
    eos::eos *eosPtr;
    eos::perfect_gas *perfectGasPtr;
    int tauCounter;

  };
  
}
#endif