TestViscousRHS.C

Go to the documentation of this file.

#include "Testing.H"
#include "Simulation.H"
#include "RKTestFixtures.H"
#include "EulerRHS.H"
#include "Stencil.H"
#include "PCPPCommUtil.H"
#include "PCPPReport.H"
#include "PCPPIntervalUtils.H"
#include "TestFixtures.H"
#include "ViscidTestFixtures.H"
#include "EulerTestFixtures.H"
#include "EulerUtil.H"
#include "ViscidUtil.H"
#include <iomanip>

using pcpp::comm::CheckResult;

typedef simulation::state::base state_t;

void TestViscidKernels(ix::test::results &serialUnitResults)
{

  typedef simulation::grid::parallel_blockstructured grid_t;

  std::cout << "TestViscidKernels.." << std::endl;
  
  //  multi-dimension test results
  std::vector<bool> dimComputeDV(2,false);
  std::vector<bool> dimComputeTV(2,false);
  std::vector<bool> dimComputeStressTensor(2,false);
  std::vector<bool> dimComputeHeatFlux(2,false);
  std::vector<bool> dimComputeDV2(2,true);
  std::vector<bool> dimComputeTEF(2,true);

  //std::vector<bool> dimUniformGridMetric(2,false);
  //std::vector<bool> dimViscidRHS(2,false);
  std::vector<bool> dimViscidUniformFlux(2,true);

  // -- Set up the grid --
  for(int nDim = 2;nDim <= 3;nDim++){
    grid_t testGrid;
    int numDim = nDim;
    std::vector<double> &realGridExtent(testGrid.PhysicalExtent());
    realGridExtent.resize(numDim*2,0.0);
    for(int iDim = 0;iDim < numDim;iDim++)
      realGridExtent[iDim*2 + 1] = 1.0;

    std::vector<size_t> &numNodes(testGrid.GridSizes());
    numNodes.resize(numDim);
    numNodes[0] = 11;
    numNodes[1] = 101;
    if(numDim == 3)
      numNodes[2] = 1001;
    
    //   numNodes[0] = 4;
    //   numNodes[1] = 4;
    //   numNodes[2] = 4;
    
    size_t numNodesDomain = 1;
    size_t numCellsDomain = 1;
    std::vector<double> dX(numDim,0.0);
    std::vector<double> gridMetric(numDim,0.0);
    std::vector<double> gridJacobian(1,1.0);
    int gridType = 0; // choose uniform rectangular metric

    for(int iDim = 0;iDim < numDim;iDim++){
      numNodesDomain *= numNodes[iDim];
      numCellsDomain *= (numNodes[iDim]-1);
      dX[iDim] = (realGridExtent[iDim*2+1] - realGridExtent[iDim*2])/(numNodes[iDim]-1);
      gridMetric[iDim] = 1.0/dX[iDim];
      gridJacobian[0] *= gridMetric[iDim];
    }
    testGrid.SetGridSpacings(dX);
    if(testGrid.Finalize()){
      std::cout << "Grid finalization failed." << std::endl;
      return;
    }

    std::vector<size_t> bufferInterval;
    testGrid.PartitionBufferInterval().Flatten(bufferInterval);

    double a = 1.0;
    double t = 1.0;
    double gamma = 1.4;
    double specGasConst = 1.0;
    double dt = 1.0;
    double cfl = 1.0;
    std::vector<double> vVec(numDim,0.0);
    vVec[0] = 1.0;
    vVec[1] = -1.0;
    if(numDim == 3) vVec[2] = 2.0;
    double v = vVec[0];
    
    // in non-dimensional, ideal gas, P=rho*R*T (R=1)
    double pressure = 1.0;
    double density = 1.2;
    double temperature = pressure/density/specGasConst;

    std::vector<double> expectedRhoM1(numNodesDomain,1.0/density);
    std::vector<double> expectedFlux((numDim+2)*numNodesDomain,0.0);
    std::vector<double> expectedPressure(numNodesDomain,pressure);
    std::vector<double> expectedTemperature(numNodesDomain,temperature);
    std::vector<double> expectedVHat(numDim*numNodesDomain,0.0);
    std::vector<double> expectedVelocity(numDim*numNodesDomain,0.0);
    std::vector<double> expectedInternalEnergy(numNodesDomain,pressure/(gamma-1.0));
    std::vector<double> expectedRHS((numDim+2)*numNodesDomain,0.0);
    
    std::vector<double> ke(numNodesDomain,0.0);
    std::vector<double> rho(numNodesDomain,density);
    // calculated below
    std::vector<double> rhoV(numDim*numNodesDomain,0.0);
    std::vector<double> rhoE(numNodesDomain,0.0); 
    // dX       = <.1,.01,.001>
    // Velocity = <1.0,-1.0,0>
    // rhoV     = 1.0 * Velocity
    // vHat     = Velocity/dX
    for(int iDim = 0;iDim < numDim;iDim++){
      for(int iPoint = 0;iPoint < numNodesDomain;iPoint++){
        rhoV[iDim*numNodesDomain + iPoint] = rho[iPoint]*vVec[iDim];
        ke[iPoint] += .5*rho[iPoint]*vVec[iDim]*vVec[iDim];
        expectedVelocity[iDim*numNodesDomain+iPoint] = vVec[iDim];
        expectedVHat[iDim*numNodesDomain+iPoint] = vVec[iDim]/dX[iDim];
      }
    } 
    // rhoE = [1.0/(gamma-1.0) + .5*rhoV*Velocity]
    for(int iPoint = 0;iPoint < numNodesDomain;iPoint++){
      //rhoE[iPoint] = 1.0/(gamma-1.0) + ke[iPoint];
      rhoE[iPoint] = expectedPressure[iPoint]/(gamma-1.0)+ ke[iPoint];
    }

    std::vector<double> T(numNodesDomain,0.0);
    std::vector<double> p(numNodesDomain,0.0);
    std::vector<double> rhom1(numNodesDomain,0.0);
    std::vector<double> velocity(numDim*numNodesDomain,-.1);
    std::vector<double> internalEnergy(numNodesDomain,0.);
    
    std::vector<double *> vBuffers(numDim,(double *)NULL);
    for(int iDim = 0;iDim < numDim;iDim++)
      vBuffers[iDim] = &velocity[iDim*numNodesDomain];

    bool computeDVWorks = true;

    T.resize(0);
    p.resize(0);
    velocity.resize(0);
    rhom1.resize(0);
    T.resize(numNodesDomain,0.0);
    p.resize(numNodesDomain,0.0);
    velocity.resize(numDim*numNodesDomain,-.1);
    rhom1.resize(numNodesDomain,0.0);
    
    std::vector<double *> myStateBuffers(numDim+2,NULL);
    myStateBuffers[0] = &rho[0];
    for(int iDim = 0;iDim < numDim;iDim++)
      myStateBuffers[iDim+1] = &rhoV[iDim*numNodesDomain];
    myStateBuffers[numDim+1] = &rhoE[0];

    std::vector<double *> myDVBuffers(4+numDim,NULL);
    myDVBuffers[0]    = &p[0];
    myDVBuffers[1]    = &T[0];
    myDVBuffers[2]    = &rhom1[0];

    for(int iDim = 0;iDim < numDim;iDim++){
      vBuffers[iDim] = &velocity[iDim*numNodesDomain];
      myDVBuffers[3+iDim] = vBuffers[iDim];
    }
    myDVBuffers[numDim+3] = &internalEnergy[0];

    pcpp::IndexIntervalType bufferExtent;
    bufferExtent.InitSimple(numNodes);

    // setup a new equation of state object for this grid
    eos::perfect_gas myEOS;
    //myEOS.InitializeBufferExtents(bufferExtent,numNodes);
    myEOS.SetupPressureBuffer(&p[0]);
    myEOS.SetupTemperatureBuffer(&T[0]);
    myEOS.SetupSpecificVolumeBuffer(&rhom1[0]);
    myEOS.SetupInternalEnergyBuffer(&internalEnergy[0]);

    myEOS.SetSpecificGasConstant(specGasConst);
    myEOS.SetGamma(gamma);

    if(myEOS.InitializeMaterialProperties()){
      std::cout << "Failed to initialize EOS" << std::endl;
    }

    if(euler::util::ComputeDVBuffer(bufferExtent,numNodes,myStateBuffers,myDVBuffers))
    {
      std::cout << "ComputeDVBuffer failed." << std::endl;
    }

    // now calculate the pressure and temperature
    if(myEOS.ComputePressureBuffer(bufferExtent,numNodes)) {
      std::cout << "ComputePressure failed." << std::endl;
    }

    if(myEOS.ComputeTemperatureBuffer(bufferExtent,numNodes)) {
      std::cout << "ComputeTemperature failed." << std::endl;
    }
    
    bool printPressure=true;
    bool printTemperature=true;
    bool printVelocity=true;
    bool printVolume=true;
    for(int iNode = 0;iNode < numNodesDomain;iNode++){
      if(std::abs(p[iNode] - expectedPressure[iNode])    > 1e-15){
        computeDVWorks = false;
        if(printPressure) {
          std::cout << "In ComputeDVTest, Expected pressure=" << expectedPressure[iNode]
                    << " got pressure=" << p[iNode] << " at iNode=" << iNode << std::endl;
          std::cout << "No further error messages will be generated for pressure" << std::endl;
          printPressure=false;
        }
      }
      if(std::abs(T[iNode] - expectedTemperature[iNode]) > 1e-15){
        computeDVWorks = false;
        if(printTemperature) {
          std::cout << "In ComputeDVTest, Expected temperature=" << expectedTemperature[iNode]
                    << " got temperature=" << T[iNode] << " at iNode=" << iNode << std::endl;
          std::cout << "No further error messages will be generated for temperature" << std::endl;
          printTemperature=false;
        }
      }
      if(std::abs(rhom1[iNode] - expectedRhoM1[iNode])   > 1e-15){
        computeDVWorks = false;
        if(printVolume) {
          std::cout << "In ComputeDVTest, Expected volume=" << expectedRhoM1[iNode]
                    << " got volume=" << rhom1[iNode] << " at iNode=" << iNode << std::endl;
          std::cout << "No further error messages will be generated for volume" << std::endl;
          printVolume=false;
        }
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        if(std::abs(velocity[iDim*numNodesDomain+iNode] - 
                    expectedVelocity[iDim*numNodesDomain+iNode]) > 1e-15){
          computeDVWorks = false;
          if(printVelocity) {
            std::cout << "In ComputeDVTest, Expected velocity=" << expectedVelocity[iNode]
                      << " got velocity=" << velocity[iNode] << " at iNode=" 
                      << iNode << " iDim=" << iDim << std::endl;
            std::cout << "No further error messages will be generated for velocity" << std::endl;
            printVelocity=false;
          }
        }
      }
    }

    //
    // test calculation of the transport coefficients
    //
    std::cout << "Running tests for "<<numDim<<"D ComputeTV" << std::endl;

    std::vector<double> mu(numNodesDomain,0.0);
    std::vector<double> lambda(numNodesDomain,0.0);
    std::vector<double> kappa(numNodesDomain,0.0);

    std::vector<double *> myTVBuffers(3,NULL);
    myTVBuffers[0]    = &mu[0];
    myTVBuffers[1]    = &lambda[0];
    myTVBuffers[2]    = &kappa[0];


    double power = 0.666;
    double beta = 5;
    double bulkViscFac = 2;
    double Cp = gamma/(gamma-1);
    double Pr = 0.75;
    double tempMu = beta*pow(expectedTemperature[0],power);
    double tempLambda = (bulkViscFac - 2.0/3.0)*tempMu;
    double tempKappa = Cp*tempMu/Pr;
    std::vector<double> expectedMu(numNodesDomain,tempMu);
    std::vector<double> expectedLambda(numNodesDomain,tempLambda);
    std::vector<double> expectedKappa(numNodesDomain,tempKappa);
    bool computeTVWorks=true;
    bool muFailed = false;
    bool lambdaFailed = false;
    bool kappaFailed = false;

    if(viscid::util::ComputeTVBufferPower(
      bufferExtent,numNodes,myDVBuffers[1],myTVBuffers,beta,power,bulkViscFac,Cp,Pr))
    {
      std::cout << "ComputeTVBuffer failed." << std::endl;
    }
    
    bool printMu = true;
    bool printLambda = true;
    bool printKappa = true;
    double error = 1;
    for(int iNode = 0;iNode < numNodesDomain;iNode++){
      error = std::abs((mu[iNode] - expectedMu[iNode])/expectedMu[iNode]);
      if(error > 1e-15){
        computeTVWorks = false;
        muFailed = true;
        if(printMu) {
          std::cout << "In ComputeTvTest, Expected mu=" << expectedMu[iNode]
                    << " got mu=" << mu[iNode] << " at iNode=" << iNode << std::endl;
          std::cout << "Temperature was " << T[iNode] << std::endl;
          std::cout << "No further error messages will be generated for mu" << std::endl;
          printMu=false;
        }
      }
      error = std::abs((lambda[iNode] - expectedLambda[iNode])/expectedLambda[iNode]);
      if(error > 1e-15){
        computeTVWorks = false;
        lambdaFailed = true;
        if(printLambda) {
          std::cout << "In ComputeTvTest, Expected lambda=" << expectedLambda[iNode]
                    << " got lambda=" << lambda[iNode] << " at iNode=" << iNode << std::endl;
          std::cout << "Temperature was " << T[iNode] << std::endl;
          std::cout << "No further error messages will be generated for lambda" << std::endl;
          printLambda=false;
        }
      }
      error = std::abs((kappa[iNode] - expectedKappa[iNode])/expectedKappa[iNode]);
      if(error > 1e-15){
        computeTVWorks = false;
        kappaFailed = true;
        if(printKappa) {
          std::cout << "In ComputeTvTest, Expected kappa=" << expectedKappa[iNode]
                    << " got kappa=" << kappa[iNode] << " at iNode=" << iNode << std::endl;
          std::cout << "Temperature was " << T[iNode] << std::endl;
          std::cout << "No further error messages will be generated for kappa" << std::endl;
          printKappa=false;
        }
      }
    }
    
    dimComputeTV[numDim-2] = computeTVWorks;
    if(muFailed){
      std::cout << "Mu failed" << std::endl;
    }
    if(lambdaFailed){
      std::cout << "Lambda failed" << std::endl;
    }
    if(kappaFailed){
      std::cout << "Kappa failed" << std::endl;
    }

    std::cout << "Done running tests for "<<numDim<<"D ComputeTV" << std::endl;
    std::cout << "Running tests for "<<numDim<<"D ComputeStressTensor" << std::endl;

    // build a velocity gradient so we can test the calculation of the stress tensor
    // in 3D...
    // | du/dx dv/dx dw/dx |   | [1..1] [2..2] [3..3]
    // | du/dy dv/dy dw/dy | = | [4..4] [5..5] [6..6]
    // | du/dz dv/dz dw/dz |   | [7..7] [8..8] [9..9]
    std::vector<double *> myVelGradBuffer(numDim,NULL);
    for(int i=0;i<numDim;i++){
      myVelGradBuffer[i] = new double [numDim*numNodesDomain];
      for(int j=0;j<numDim;j++){
        int offset=j*numNodesDomain;
        double value = static_cast<double>(i*numDim+j+1);
        for(int k=0;k<numNodesDomain;k++){
          //if(k==0){
            //std::cout << "i, j, velGrad" << i << " " << j << " " << value << std::endl;
          //}
          myVelGradBuffer[i][offset+k]=value;
        }
      }
    }

    // calculate the stress tensor
    // the stress tensor is symmetric and stored in memory like (in 3D)
    // [ {tau11} {tau12} {tau13} {tau22} {tau23} {tau33} ]
    // myTauBuffer has pointers into memory for all numDim^2 components
    // [ {tau11} {tau12} {tau13} {tau22} {tau23} {tau33} ]
    std::vector<double *> myTauBuffer(numDim*numDim,NULL);
    std::vector<double> myTau(numDim*(numDim+1)*numNodesDomain/2,0);
    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){
          myTauBuffer[index] = &myTau[bufferIndex];
          bufferIndex += numNodesDomain;
        } else {
          myTauBuffer[index] = myTauBuffer[j*numDim+i];
        }
      }
    }

    double Re = 100.;
    if(viscid::util::ComputeTauBuffer(bufferExtent, numNodes, gridType, gridMetric, gridJacobian, Re, 
                                      myTVBuffers, myVelGradBuffer, myTauBuffer)){
      std::cout << "Calculation of the stress tensor failed" << std::endl;
    }

    // Construct the expected stress tensor
    // tau_ij = mu/Re * (dv_i/dx_j+dv_j/dx_i) + lambda/Re * (dv_k/dx_k) delta_ij
    // Re = 100, mu and lambda set above from T=3.5
    // gridmetrics are applied by ComputeStressTensor, so account for this here
    std::vector<double> tauExpected;

    if(numDim == 2){
      tauExpected.push_back(tempMu/Re*(1*gridMetric[0]+1*gridMetric[0]) 
        + tempLambda/Re*(1*gridMetric[0]+4*gridMetric[1]));                 //tau_11
      tauExpected.push_back(tempMu/Re*(2*gridMetric[0]+3*gridMetric[1]));   //tau_12
      tauExpected.push_back(tempMu/Re*(2*gridMetric[0]+3*gridMetric[1]));   //tau_21 = tau_12
      tauExpected.push_back(tempMu/Re*(4*gridMetric[1]+4*gridMetric[1]) 
        + tempLambda/Re*(1*gridMetric[0]+4*gridMetric[1]));                 //tau_22

      

    } else {

    // | du/dxi   dv/dxi   dw/dxi   |   | [1..1] [2..2] [3..3]
    // | du/deta  dv/deta  dw/deta  | = | [4..4] [5..5] [6..6]
    // | du/dzeta dv/dzeta dw/dzeta |   | [7..7] [8..8] [9..9]
    // dxidx   = metric[0]
    // detady  = metric[1]
    // dzetadz = metric[2]

      tauExpected.push_back(tempMu/Re*(1*gridMetric[0]+1*gridMetric[0]) 
                            + tempLambda/Re*(1*gridMetric[0]+5*gridMetric[1]+9*gridMetric[2])); //tau_11
      tauExpected.push_back(tempMu/Re*(2*gridMetric[0]+4*gridMetric[1]));   //tau_12
      tauExpected.push_back(tempMu/Re*(3*gridMetric[0]+7*gridMetric[2]));   //tau_13
      tauExpected.push_back(tempMu/Re*(2*gridMetric[0]+4*gridMetric[1]));   //tau_21 = tau12
      tauExpected.push_back(tempMu/Re*(5*gridMetric[1]+5*gridMetric[1]) 
                            + tempLambda/Re*(1*gridMetric[0]+5*gridMetric[1]+9*gridMetric[2])); //tau_22
      tauExpected.push_back(tempMu/Re*(6*gridMetric[1]+8*gridMetric[2]));   //tau_23
      tauExpected.push_back(tempMu/Re*(3*gridMetric[0]+7*gridMetric[2]));   //tau_31 = tau13
      tauExpected.push_back(tempMu/Re*(6*gridMetric[1]+8*gridMetric[2]));   //tau_32 = tau23
      tauExpected.push_back(tempMu/Re*(9*gridMetric[2]+9*gridMetric[2]) 
                            + tempLambda/Re*(1*gridMetric[0]+5*gridMetric[1]+9*gridMetric[2])); //tau_33
    }

    for(int i=0;i<tauExpected.size();i++){
      tauExpected[i] *= gridJacobian[0];
    }

    // Compare expected against actual
    bool tauFailed=false;
    bool computeTauWorks=true;
    for(int i=0;i<numDim*numDim;i++){
      bool printTau=true;
      for(int iNode = 0;iNode < numNodesDomain;iNode++){
        double error = std::abs((myTauBuffer[i][iNode] - tauExpected[i])/tauExpected[i]);
        if(error > 1e-15){
          computeTauWorks = false;
          tauFailed = true;
          if(printTau) {
            std::cout << "In ComputeTauTest, Expected tau[" << i << "]=" << tauExpected[i]
                      << " got tau=" << myTauBuffer[i][iNode] << " at iNode=" << iNode 
                      << " error=" << error << std::endl;
            std::cout << "No further error messages will be generated for this tau index" 
                      << std::endl;
            printTau=false;
          }
        }
      }
    }
    dimComputeStressTensor[numDim-2] = computeTauWorks;
    if(tauFailed){
      std::cout << "Tau failed" << std::endl;
    }
    std::cout << "Done running tests for "<<numDim<<"D ComputeStressTensor" << std::endl;

    std::cout << "Running tests for "<<numDim<<"D ComputeHeatFlux" << std::endl;

    // build a temperature gradient so we can test the calculation of the heat flux vector
    // in 3D...
    // | dT/dx dT/dy dT/dz |   | [1..1] [2..2] [3..3]
    std::vector<double *> myTempGradBuffer(numDim,NULL);
    std::vector<double> myTempGrad(numDim*numNodesDomain,0);
    for(int i=0;i<numDim;i++){
      myTempGradBuffer[i] = &myTempGrad[i*numNodesDomain];
      double value = static_cast<double>(i+1);
      for(int k=0;k<numNodesDomain;k++){
        myTempGradBuffer[i][k]=value;
      }
    }

    // calculate the heat flux vector
    std::vector<double *> myHeatFluxBuffer(numDim,NULL);
    std::vector<double> myHeatFlux(numDim*numNodesDomain,0);
    for(int i=0;i<numDim;i++){
      myHeatFluxBuffer[i] = &myHeatFlux[i*numNodesDomain];
    }

    if(viscid::util::ComputeHeatFluxBuffer(bufferExtent, numNodes, gridType,gridMetric, 
                                           gridJacobian, Re, myTVBuffers, myTempGradBuffer, 
                                           myHeatFluxBuffer)){
      std::cout << "Calculation of the heat flux vector failed" << std::endl;
    }

    // Construct the expected heat flux vector
    // q_i = mu/Re/Pr * (dT/dx_i) 
    // Re = 100, Pr = 0.75, mu and lambda set above from T=3.5
    // gridmetrics are applied by ComputeHeatFlux, so account for this here
    std::vector<double> heatFluxExpected;
    for(int iDim = 0;iDim < numDim;iDim++)
      heatFluxExpected.push_back(Cp*tempMu/Pr/Re*(iDim+1)*gridMetric[iDim]*gridJacobian[0]);

    // Compare expected against actual
    bool heatFluxFailed=false;
    bool computeHeatFluxWorks=true;
    for(int i=0;i<numDim;i++){
      bool printHeatFlux=true;
      for(int iNode = 0;iNode < numNodesDomain;iNode++){
        double error = std::abs((myHeatFluxBuffer[i][iNode] - heatFluxExpected[i])/
          heatFluxExpected[i]);
        if(error > 1e-15){
          computeHeatFluxWorks = false;
          heatFluxFailed = true;
          if(printHeatFlux) {
            std::cout << "In ComputeHeatFluxTest, Expected heatFlux[" << i << "]=" 
                      << heatFluxExpected[i]
                      << " got heatFlux=" << myHeatFluxBuffer[i][iNode] << " at iNode=" 
                      << iNode << " error=" << error << std::endl;
            std::cout << "No further error messages will be generated for this q index" 
                      << std::endl;
            printHeatFlux=false;
          }
        }
      }
    }
    dimComputeHeatFlux[numDim-2] = computeHeatFluxWorks;
    if(heatFluxFailed){
      std::cout << "Heat flux failed" << std::endl;
    }
    std::cout << "Done running tests for " << numDim 
              << "D ComputeHeatFlux" << std::endl;

    //
    // now test the computation of the viscous flux terms
    // This is completely fabricated and non-physical as the velocity
    // is inconsistent with the velocity gradient.
    //
    std::cout << "Running tests for " << numDim 
              << "D viscous flux" << std::endl;
    
    std::vector<size_t> fortranBufferInterval(2*numDim,1);
    for(int iDim = 0;iDim < numDim;iDim++){
      fortranBufferInterval[2*iDim+1] = bufferInterval[2*iDim+1] + 1;
    }
    
    for(int iDim = 0;iDim < numDim;iDim++){
      vBuffers[iDim] = &velocity[iDim*numNodesDomain];
    }
    
    // New sequence:
    // Compute v_i, tau_i_j and -q_i (done above in original test)
    // Compute the energy flux vector Q_i = (v_j * tau_i_j + q_i)  
    // Use new flux routine (viscidstrongflux1d) to get viscous fluxes:
    // momentum terms = metric * tau
    // energy term = metric * Q
    
    for(int iDim = 0;iDim < numDim;iDim++){
      // Compute q_i += (v_j * tau_i_j)
      for(int iVel = 0;iVel < numDim;iVel++){
        int tensorIndex = iDim*numDim + iVel;
        size_t dirOffset = iVel*numNodesDomain;

        FC_MODULE(operators,ywxpy,OPERATORS,YWXPY)
          (&numDim,&numNodesDomain,&numNodes[0],&fortranBufferInterval[0],
           &velocity[dirOffset],myTauBuffer[tensorIndex],myHeatFluxBuffer[iDim]);
      }
    }

    // Check expected v * tau
    for(int iDim = 0;iDim < numDim;iDim++){
      bool triggered = false;
      double expectedValue = 0;
      for(int iVel = 0;iVel < numDim;iVel++){
        int tensorIndex = iDim*numDim + iVel;
        expectedValue += (vVec[iVel] * tauExpected[tensorIndex]);
      }
      expectedValue += heatFluxExpected[iDim];
      for(size_t iPoint = 0;iPoint < numNodesDomain;iPoint++){
        if((std::abs(myHeatFluxBuffer[iDim][iPoint] - expectedValue)/std::abs(expectedValue)) >
           1e-12){
          dimComputeTEF[nDim-2] = false;
          if(!triggered){
            std::cout << "Failed in " << (iDim+1) << "/" << numDim 
                      << " dimension " << std::endl;
            std::cout << "TEF(" << iPoint << ") [r,e] = [" << myHeatFluxBuffer[iDim][iPoint]
                      << "," << expectedValue << "]" << std::endl;
            triggered = true;
          }
        }
      }
    }
    
    for(int iDim = 0;iDim < numDim;iDim++){

      std::vector<double> viscidFlux(numNodesDomain*(numDim+2),0.0);
      std::vector<double> expectedFlux((numDim+2)*numNodesDomain,0.0);
      std::vector<double> scaledTau(numNodesDomain,0.0);

      for(size_t iPoint = 0;iPoint < numNodesDomain;iPoint++){
        size_t fOffset = 0;

        for(int ivDim = 0;ivDim < numDim;ivDim++){
          size_t iOffset = ivDim*numNodesDomain;
          fOffset = iOffset + numNodesDomain;
          size_t vectorOffset = ivDim*numNodesDomain;

          // momentum
          int tauIndex=iDim*numDim+ivDim;
          expectedFlux[iPoint+fOffset] += tauExpected[tauIndex]*gridMetric[iDim];

          // energy
          expectedFlux[iPoint+(numDim+1)*numNodesDomain] += 
            tauExpected[tauIndex]*velocity[vectorOffset+iPoint]*gridMetric[iDim];
        }

        // thermal conduction
        expectedFlux[iPoint+(numDim+1)*numNodesDomain] +=
          heatFluxExpected[iDim]*gridMetric[iDim];

        fOffset += numNodesDomain;
      }

      int fluxDim = iDim + 1;
      int tau1=iDim*numDim;
      int tau2=tau1+1;
      int tau3=tau2+1;
      
      FC_MODULE(viscid,strongflux1d,VISCID,STRONGFLUX1D)
        (&numDim,&fluxDim,&numNodes[0],&numNodesDomain,&fortranBufferInterval[0],
         &gridType,&gridMetric[0],&myTau[0],&myHeatFlux[0],&viscidFlux[0]);

      std::vector<bool> uniformFlux(numDim+2,true);
      for(int fDim = 0 ;fDim < (numDim+2);fDim++){
        size_t fOffset = fDim * numNodesDomain;
        bool triggered = false;
        for(size_t iPoint = 0;iPoint < numNodesDomain;iPoint++){
          double error= std::abs((viscidFlux[fOffset+iPoint] - expectedFlux[fOffset+iPoint])/
            expectedFlux[fOffset+iPoint]);
          if(error > 1e-10){
            if(!triggered ){
              std::cout << "iDim " << iDim << " Flux Dimension " << fDim << std::endl;
              std::cout << std::scientific << std::setprecision(16)
                        << "\tActual:     " << viscidFlux[fOffset+iPoint] << std::endl
                        << "\tExpected:   " << expectedFlux[fOffset+iPoint] << std::endl
                        << "\tError:      " << error << std::endl; 
              triggered = true;
            }
            uniformFlux[fDim] = false;
          }
        }
      }
      for(int fDim = 0;fDim < (numDim+2);fDim++){
        if(!uniformFlux[fDim]){
          dimViscidUniformFlux[numDim-2] = false;
          std::cout << "uniform flux failed in dim,fdim " << iDim 
                    << "," << fDim << std::endl;
        }
      }
      std::cout << "Done running tests for "<<numDim<<"D viscous flux" << std::endl;

    } // loop over iDim
  }
    
  bool computeDVWorks = true;
  bool computeTVWorks = true;
  bool computeStressTensorWorks = true;
  bool computeHeatFluxWorks = true;
  bool viscidUniformFluxWorks = true;
  bool computeTotalEnergyFlux = true;
  // Make sure these tests worked in 2d and 3d
  for(int nDim = 2;nDim <= 3;nDim++){
    // compute DV tests are temporay until I have computeTV working
    computeDVWorks = computeDVWorks && dimComputeDV2[nDim-2];
    computeTVWorks = computeTVWorks && dimComputeTV[nDim-2];
    computeTotalEnergyFlux = computeTotalEnergyFlux && dimComputeTEF[nDim-2];
    computeStressTensorWorks = computeStressTensorWorks && dimComputeStressTensor[nDim-2];
    computeHeatFluxWorks = computeHeatFluxWorks && dimComputeHeatFlux[nDim-2];
    viscidUniformFluxWorks = viscidUniformFluxWorks && dimViscidUniformFlux[nDim-2];
  }
  
  serialUnitResults.UpdateResult("Viscid:Kernels:ComputeTEF",computeTotalEnergyFlux);
  serialUnitResults.UpdateResult("Viscid:Kernels:DV",computeDVWorks);
  serialUnitResults.UpdateResult("Viscid:Kernels:TV",computeTVWorks);
  serialUnitResults.UpdateResult("Viscid:Kernels:StressTensor",computeStressTensorWorks);
  serialUnitResults.UpdateResult("Viscid:Kernels:HeatFlux",computeHeatFluxWorks);
  serialUnitResults.UpdateResult("Viscid:Kernels:UniformViscidFlux",viscidUniformFluxWorks);

};
void TestViscidKernelsMetrics(ix::test::results &serialUnitResults)
{

  typedef simulation::grid::parallel_blockstructured grid_t;
  const std::string functionName("TestViscidKernelsMetrics");
  
  std::cout << "Running " << functionName << std::endl;

  bool testResult = true;

  //  multi-dimension test results
  std::vector<bool> dimComputeDV(2,false);
  std::vector<bool> dimComputeTV(2,false);
  std::vector<bool> dimComputeStressTensor(2,false);
  std::vector<bool> dimComputeHeatFlux(2,false);
  std::vector<bool> dimComputeDV2(2,true);
  std::vector<bool> dimComputeTEF(2,true);
  
  std::vector<bool> dimViscidUniformFlux(2,true);
  
  std::vector<bool> metricSuccess(simulation::grid::NUMGRIDTYPES,true);

  for(int gridType = simulation::grid::UNIRECT;
      gridType < simulation::grid::NUMGRIDTYPES;
      gridType++){
    std::cout << "Testing viscid kernels with "
              << (gridType <= simulation::grid::UNIRECT ? " Uniform-rectangular " :
                  gridType == simulation::grid::RECTILINEAR ? " Rectilinear " :
                  " Curvilinear ") << " metrics." << std::endl;

    // -- Set up the grid --
    for(int nDim = 2;nDim <= 3;nDim++){


      std::cout << "Testing viscid kernels in " << nDim << " dimensions." << std::endl;
      grid_t testGrid;
      int numDim = nDim;

      std::vector<double> &realGridExtent(testGrid.PhysicalExtent());
      realGridExtent.resize(numDim*2,0.0);
      for(int iDim = 0;iDim < numDim;iDim++)
        realGridExtent[iDim*2 + 1] = 1.0;
      
      std::vector<size_t> &numNodes(testGrid.GridSizes());
      numNodes.resize(numDim);
      numNodes[0] = 11;
      numNodes[1] = 101;
      if(numDim == 3)
        numNodes[2] = 1001;
      
      size_t numNodesDomain = 1;
      size_t numCellsDomain = 1;
      for(int iDim = 0;iDim < numDim;iDim++){
        numNodesDomain *= numNodes[iDim];
        numCellsDomain *= (numNodes[iDim]-1);
      }
      
      std::vector<double> dX;
      std::vector<double> gridMetric;
      std::vector<double> gridJacobian;
      
      if(gridType <= simulation::grid::UNIRECT){
        dX.resize(numDim,0.0);
        gridMetric.resize(numDim,0.0);
        gridJacobian.resize(1,1.0);
        for(int iDim = 0;iDim < numDim;iDim++){
          dX[iDim] = (realGridExtent[iDim*2+1] - realGridExtent[iDim*2])/(numNodes[iDim]-1);
          gridMetric[iDim] = 1.0/dX[iDim];
          gridJacobian[0] *= gridMetric[iDim];
        }
        testGrid.SetGridSpacings(dX);
      } else if(gridType == simulation::grid::RECTILINEAR){
        dX.resize(numDim*numNodesDomain,0.0);
        gridMetric.resize(numDim*numNodesDomain,0.0);
        gridJacobian.resize(numNodesDomain,1.0);
        for(int iDim = 0;iDim < numDim;iDim++){
          size_t dimOffset = iDim*numNodesDomain;
          for(size_t iPoint = 0;iPoint < numNodesDomain;iPoint++){
            dX[dimOffset+iPoint] = (realGridExtent[iDim*2+1] - realGridExtent[iDim*2])/(numNodes[iDim]-1);
            gridMetric[dimOffset+iPoint] = 1.0/dX[dimOffset+iPoint];
            gridJacobian[iPoint] *= gridMetric[dimOffset+iPoint];
          }
        }
      } else {
        dX.resize(numDim*numDim*numNodesDomain,0.0);
        gridMetric.resize(numDim*numDim*numNodesDomain,0.0);
        gridJacobian.resize(numNodesDomain,0.0);
        for(int iDim = 0;iDim < numDim;iDim++){
          size_t iIndex = iDim*numDim;
          for(int jDim = 0;jDim < numDim;jDim++){
            size_t ijIndex = iIndex + jDim;
            size_t metricOffset = ijIndex*numNodesDomain;
            if(iDim == jDim){
              for(int iPoint = 0;iPoint < numNodesDomain;iPoint++){
                dX[metricOffset+iPoint] = (realGridExtent[iDim*2+1] - realGridExtent[iDim*2])/(numNodes[iDim]-1);
                gridMetric[metricOffset+iPoint] = 1.0/dX[metricOffset+iPoint];
                gridJacobian[iPoint] *= gridMetric[metricOffset+iPoint];
              }
            }
          }
        }
      }
      if(testGrid.Finalize()){
        std::cout << "Grid finalization failed." << std::endl;
        return;
      }
      
      std::vector<size_t> bufferInterval;
      testGrid.PartitionBufferInterval().Flatten(bufferInterval);
    
      double a = 1.0;
      double t = 1.0;
      double gamma = 1.4;
      double specGasConst = 1.0;
      double dt = 1.0;
      double cfl = 1.0;
      std::vector<double> vVec(numDim,0.0);
      vVec[0] = 1.0;
      vVec[1] = -1.0;
      if(numDim == 3) vVec[2] = 2.0;
      double v = vVec[0];
    
      // in non-dimensional, ideal gas, P=rho*R*T (R=1)
      double pressure = 1.0;
      double density = 1.2;
      double temperature = pressure/density/specGasConst;
    
      std::vector<double> expectedRhoM1(numNodesDomain,1.0/density);
      std::vector<double> expectedFlux((numDim+2)*numNodesDomain,0.0);
      std::vector<double> expectedPressure(numNodesDomain,pressure);
      std::vector<double> expectedTemperature(numNodesDomain,temperature);
      std::vector<double> expectedVHat(numDim*numNodesDomain,0.0);
      std::vector<double> expectedVelocity(numDim*numNodesDomain,0.0);
      std::vector<double> expectedRHS((numDim+2)*numNodesDomain,0.0);
    
      std::vector<double> ke(numNodesDomain,0.0);
      std::vector<double> rho(numNodesDomain,density);
      std::vector<double> rhoV(numDim*numNodesDomain,1.0);
      std::vector<double> rhoE(numNodesDomain,0.0); // calculated below
      // dX       = <.1,.01,.001>
      // Velocity = <1.0,-1.0,0>
      // rhoV     = 1.0 * Velocity
      // vHat     = V * xiHat

      for(int iDim = 0;iDim < numDim;iDim++){
        double dx = dX[iDim];
        if(gridType == simulation::grid::RECTILINEAR){
          dx = dX[iDim*numNodesDomain];
        } else if (gridType == simulation::grid::CURVILINEAR) {
          dx = dX[iDim*numDim*numNodesDomain];
        }
        for(int iPoint = 0;iPoint < numNodesDomain;iPoint++){
          rhoV[iDim*numNodesDomain + iPoint] = rho[iPoint]*vVec[iDim];
          ke[iPoint] += .5*rho[iPoint]*vVec[iDim]*vVec[iDim];
          expectedVelocity[iDim*numNodesDomain+iPoint] = vVec[iDim];
          expectedVHat[iDim*numNodesDomain+iPoint] = vVec[iDim]/dx;
        }
      } 
      // rhoE = [1.0/(gamma-1.0) + .5*rhoV*Velocity]
      for(int iPoint = 0;iPoint < numNodesDomain;iPoint++){
        //rhoE[iPoint] = 1.0/(gamma-1.0) + ke[iPoint];
        rhoE[iPoint] = expectedPressure[iPoint]/(gamma-1.0)+ ke[iPoint];
      }
    
      std::vector<double> T(numNodesDomain,0.0);
      std::vector<double> p(numNodesDomain,0.0);
      std::vector<double> rhom1(numNodesDomain,0.0);
      std::vector<double> velocity(numDim*numNodesDomain,-.1);
      std::vector<double> internalEnergy(numNodesDomain,0.);
    
      std::vector<double *> vBuffers(numDim,(double *)NULL);
      for(int iDim = 0;iDim < numDim;iDim++)
        vBuffers[iDim] = &velocity[iDim*numNodesDomain];
    
      bool computeDVWorks = true;
      bool pfailed = false;
      bool tfailed = false;
      bool rfailed = false;
      bool vfailed = false;
    
      T.resize(0);
      p.resize(0);
      velocity.resize(0);
      rhom1.resize(0);
      T.resize(numNodesDomain,0.0);
      p.resize(numNodesDomain,0.0);
      velocity.resize(numDim*numNodesDomain,-.1);
      rhom1.resize(numNodesDomain,0.0);
    
      std::vector<double *> myStateBuffers(numDim+2,NULL);
      myStateBuffers[0] = &rho[0];
      for(int iDim = 0;iDim < numDim;iDim++)
        myStateBuffers[iDim+1] = &rhoV[iDim*numNodesDomain];
      myStateBuffers[numDim+1] = &rhoE[0];
    
      std::vector<double *> myDVBuffers(4+numDim,NULL);
      myDVBuffers[0]    = &p[0];
      myDVBuffers[1]    = &T[0];
      myDVBuffers[2]    = &rhom1[0];
    
      for(int iDim = 0;iDim < numDim;iDim++){
        vBuffers[iDim] = &velocity[iDim*numNodesDomain];
        myDVBuffers[3+iDim] = vBuffers[iDim];
      }
      myDVBuffers[numDim+3] = &internalEnergy[0];
    
      pcpp::IndexIntervalType bufferExtent;
      bufferExtent.InitSimple(numNodes);

      // setup a new equation of state object for this grid
      eos::perfect_gas myEOS;
      //myEOS.InitializeBufferExtents(bufferExtent,numNodes);
      myEOS.SetupPressureBuffer(&p[0]);
      myEOS.SetupTemperatureBuffer(&T[0]);
      myEOS.SetupSpecificVolumeBuffer(&rhom1[0]);
      myEOS.SetupInternalEnergyBuffer(&internalEnergy[0]);
  
      myEOS.SetSpecificGasConstant(specGasConst);
      myEOS.SetGamma(gamma);
  
      if(myEOS.InitializeMaterialProperties()){
        std::cout << "Failed to initialize EOS" << std::endl;
      }
    
      if(euler::util::ComputeDVBuffer(bufferExtent,numNodes,myStateBuffers,myDVBuffers))
      {
        std::cout << "ComputeDVBuffer failed." << std::endl;
      }

      // now calculate the pressure and temperature
      if(myEOS.ComputePressureBuffer(bufferExtent,numNodes)) {
        std::cout << "ComputePressure failed." << std::endl;
      }
  
      if(myEOS.ComputeTemperatureBuffer(bufferExtent,numNodes)) {
        std::cout << "ComputeTemperature failed." << std::endl;
      }
    

      bool printPressure=true;
      bool printTemperature=true;
      bool printVelocity=true;
      bool printVolume=true;
      for(int iNode = 0;iNode < numNodesDomain;iNode++){
        if(std::abs(p[iNode] - expectedPressure[iNode])    > 1e-15){
          computeDVWorks = false;
          if(printPressure) {
            std::cout << "In ComputeDVTest, Expected pressure=" << expectedPressure[iNode]
                      << " got pressure=" << p[iNode] << " at iNode=" << iNode << std::endl;
            std::cout << "No further error messages will be generated for pressure" << std::endl;
            printPressure=false;
          }
        }
        if(std::abs(T[iNode] - expectedTemperature[iNode]) > 1e-15){
          computeDVWorks = false;
          if(printTemperature) {
            std::cout << "In ComputeDVTest, Expected temperature=" << expectedTemperature[iNode]
                      << " got temperature=" << T[iNode] << " at iNode=" << iNode << std::endl;
            std::cout << "No further error messages will be generated for temperature" << std::endl;
            printTemperature=false;
          }
        }
        if(std::abs(rhom1[iNode] - expectedRhoM1[iNode])   > 1e-15){
          computeDVWorks = false;
          if(printVolume) {
            std::cout << "In ComputeDVTest, Expected volume=" << expectedRhoM1[iNode]
                      << " got volume=" << rhom1[iNode] << " at iNode=" << iNode << std::endl;
            std::cout << "No further error messages will be generated for volume" << std::endl;
            printVolume=false;
          }
        }
        for(int iDim = 0;iDim < numDim;iDim++){
          if(std::abs(velocity[iDim*numNodesDomain+iNode] - 
                      expectedVelocity[iDim*numNodesDomain+iNode]) > 1e-15){
            computeDVWorks = false;
            if(printVelocity) {
              std::cout << "In ComputeDVTest, Expected velocity=" << expectedVelocity[iNode]
                        << " got velocity=" << velocity[iNode] << " at iNode=" 
                        << iNode << " iDim=" << iDim << std::endl;
              std::cout << "No further error messages will be generated for velocity" << std::endl;
              printVelocity=false;
            }
          }
        }
      }
    
      dimComputeDV2[numDim-2] = computeDVWorks ;

      //
      // test calculation of the transport coefficients
      //
      std::cout << "Running tests for "<<numDim<<"D ComputeTV" << std::endl;

      std::vector<double> mu(numNodesDomain,0.0);
      std::vector<double> lambda(numNodesDomain,0.0);
      std::vector<double> kappa(numNodesDomain,0.0);

      std::vector<double *> myTVBuffers(3,NULL);
      myTVBuffers[0]    = &mu[0];
      myTVBuffers[1]    = &lambda[0];
      myTVBuffers[2]    = &kappa[0];


      double power = 0.666;
      double beta = 5;
      double bulkViscFac = 2;
      double Cp = gamma/(gamma-1);
      double Pr = 0.75;
      double tempMu = beta*pow(expectedTemperature[0],power);
      double tempLambda = (bulkViscFac - 2.0/3.0)*tempMu;
      double tempKappa = Cp*tempMu/Pr;
      std::vector<double> expectedMu(numNodesDomain,tempMu);
      std::vector<double> expectedLambda(numNodesDomain,tempLambda);
      std::vector<double> expectedKappa(numNodesDomain,tempKappa);
      bool computeTVWorks=true;
      bool muFailed = false;
      bool lambdaFailed = false;
      bool kappaFailed = false;

      if(viscid::util::ComputeTVBufferPower(
                                            bufferExtent,numNodes,myDVBuffers[1],myTVBuffers,beta,power,bulkViscFac,Cp,Pr))
        {
          std::cout << "ComputeTVBuffer failed." << std::endl;
        }
    
      bool printMu = true;
      bool printLambda = true;
      bool printKappa = true;
      double error = 1;
      for(int iNode = 0;iNode < numNodesDomain;iNode++){
        error = std::abs((mu[iNode] - expectedMu[iNode])/expectedMu[iNode]);
        if(error > 1e-15){
          computeTVWorks = false;
          muFailed = true;
          if(printMu) {
            std::cout << "In ComputeTvTest, Expected mu=" << expectedMu[iNode]
                      << " got mu=" << mu[iNode] << " at iNode=" << iNode << std::endl;
            std::cout << "Temperature was " << T[iNode] << std::endl;
            std::cout << "No further error messages will be generated for mu" << std::endl;
            printMu=false;
          }
        }
        error = std::abs((lambda[iNode] - expectedLambda[iNode])/expectedLambda[iNode]);
        if(error > 1e-15){
          computeTVWorks = false;
          lambdaFailed = true;
          if(printLambda) {
            std::cout << "In ComputeTvTest, Expected lambda=" << expectedLambda[iNode]
                      << " got lambda=" << lambda[iNode] << " at iNode=" << iNode << std::endl;
            std::cout << "Temperature was " << T[iNode] << std::endl;
            std::cout << "No further error messages will be generated for lambda" << std::endl;
            printLambda=false;
          }
        }
        error = std::abs((kappa[iNode] - expectedKappa[iNode])/expectedKappa[iNode]);
        if(error > 1e-15){
          computeTVWorks = false;
          kappaFailed = true;
          if(printKappa) {
            std::cout << "In ComputeTvTest, Expected kappa=" << expectedKappa[iNode]
                      << " got kappa=" << kappa[iNode] << " at iNode=" << iNode << std::endl;
            std::cout << "Temperature was " << T[iNode] << std::endl;
            std::cout << "No further error messages will be generated for kappa" << std::endl;
            printKappa=false;
          }
        }
      }
    
      dimComputeTV[numDim-2] = computeTVWorks;
      if(muFailed){
        std::cout << "Mu failed" << std::endl;
      }
      if(lambdaFailed){
        std::cout << "Lambda failed" << std::endl;
      }
      if(kappaFailed){
        std::cout << "Kappa failed" << std::endl;
      }

      std::cout << "Done running tests for "<<numDim<<"D ComputeTV" << std::endl;
      std::cout << "Running tests for "<<numDim<<"D ComputeStressTensor" << std::endl;

      // build a velocity gradient so we can test the calculation of the stress tensor
      // in 3D...
      // | du/dxi   dv/dxi   dw/dxi   |   | [1..1] [2..2] [3..3]
      // | du/deta  dv/deta  dw/deta  | = | [4..4] [5..5] [6..6]
      // | du/dzeta dv/dzeta dw/dzeta |   | [7..7] [8..8] [9..9]
      std::vector<double *> myVelGradBuffer(numDim,NULL);
      for(int i=0;i<numDim;i++){
        myVelGradBuffer[i] = new double [numDim*numNodesDomain];
        for(int j=0;j<numDim;j++){
          int offset=j*numNodesDomain;
          double value = static_cast<double>(i*numDim+j+1);
          for(int k=0;k<numNodesDomain;k++){
            //if(k==0){
            //std::cout << "i, j, velGrad" << i << " " << j << " " << value << std::endl;
            //}
            myVelGradBuffer[i][offset+k]=value;
          }
        }
      }

      // calculate the stress tensor
      // the stress tensor is symmetric and stored in memory like (in 3D)
      // [ {tau11} {tau12} {tau13} {tau22} {tau23} {tau33} ]
      // myTauBuffer has pointers into memory for all numDim^2 components
      // [ {tau11} {tau12} {tau13} {tau22} {tau23} {tau33} ]
      std::vector<double *> myTauBuffer(numDim*numDim,NULL);
      std::vector<double> myTau(numDim*(numDim+1)*numNodesDomain/2,0);
      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){
            myTauBuffer[index] = &myTau[bufferIndex];
            bufferIndex += numNodesDomain;
          } else {
            myTauBuffer[index] = myTauBuffer[j*numDim+i];
          }
        }
      }

      double Re = 100.;
      if(viscid::util::ComputeTauBuffer(bufferExtent, numNodes, gridType, gridMetric, gridJacobian, Re, 
                                        myTVBuffers, myVelGradBuffer, myTauBuffer)){
        std::cout << "Calculation of the stress tensor failed" << std::endl;
      }

      // Construct the expected stress tensor
      // tau_ij = mu/Re * (dv_i/dx_j+dv_j/dx_i) + lambda/Re * (dv_k/dx_k) delta_ij
      // Re = 100, mu and lambda set above from T=3.5
      // gridmetrics are applied by ComputeStressTensor, so account for this here
      std::vector<double> tauExpected;

      if(numDim == 2){

        double metricX = gridMetric[0];
        double metricY = gridMetric[1];

        if(gridType == simulation::grid::RECTILINEAR) {
          metricY = gridMetric[numNodesDomain];
        } else if (gridType == simulation::grid::CURVILINEAR) {
          metricY = gridMetric[3*numNodesDomain];
        }
      
        tauExpected.push_back(tempMu/Re*(1*metricX+1*metricX) 
                              + tempLambda/Re*(1*metricX+4*metricY));//tau_11
        tauExpected.push_back(tempMu/Re*(2*metricX+3*metricY));      //tau_12
        tauExpected.push_back(tempMu/Re*(2*metricX+3*metricY));      //tau_21 = tau_12
        tauExpected.push_back(tempMu/Re*(4*metricY+4*metricY) 
                              + tempLambda/Re*(1*metricX+4*metricY));//tau_22

      } else {

        // | du/dxi   dv/dxi   dw/dxi   |   | [1..1] [2..2] [3..3]
        // | du/deta  dv/deta  dw/deta  | = | [4..4] [5..5] [6..6]
        // | du/dzeta dv/dzeta dw/dzeta |   | [7..7] [8..8] [9..9]
        // dxidx   = metric[0]
        // detady  = metric[1]
        // dzetadz = metric[2]
        double metricX = gridMetric[0];
        double metricY = gridMetric[1];
        double metricZ = gridMetric[2];
        if(gridType == simulation::grid::RECTILINEAR){
          metricY = gridMetric[numNodesDomain];
          metricZ = gridMetric[2*numNodesDomain];
        } else if (gridType == simulation::grid::CURVILINEAR){
          metricY = gridMetric[4*numNodesDomain];
          metricZ = gridMetric[8*numNodesDomain];
        }

        tauExpected.push_back(tempMu/Re*(1*metricX+1*metricX) 
                              + tempLambda/Re*(1*metricX+5*metricY+9*metricZ)); //tau_11
        tauExpected.push_back(tempMu/Re*(2*metricX+4*metricY));   //tau_12
        tauExpected.push_back(tempMu/Re*(3*metricX+7*metricZ));   //tau_13
        tauExpected.push_back(tempMu/Re*(2*metricX+4*metricY));   //tau_21 = tau12
        tauExpected.push_back(tempMu/Re*(5*metricY+5*metricY) 
                              + tempLambda/Re*(1*metricX+5*metricY+9*metricZ)); //tau_22
        tauExpected.push_back(tempMu/Re*(6*metricY+8*metricZ));   //tau_23
        tauExpected.push_back(tempMu/Re*(3*metricX+7*metricZ));   //tau_31 = tau13
        tauExpected.push_back(tempMu/Re*(6*metricY+8*metricZ));   //tau_32 = tau23
        tauExpected.push_back(tempMu/Re*(9*metricZ+9*metricZ) 
                              + tempLambda/Re*(1*metricX+5*metricY+9*metricZ)); //tau_33
      }

      for(int i=0;i<tauExpected.size();i++){
        tauExpected[i] *= gridJacobian[0];
      }

      // Compare expected against actual
      bool tauFailed=false;
      bool computeTauWorks=true;
      for(int i=0;i<numDim*numDim;i++){
        bool printTau=true;
        for(int iNode = 0;iNode < numNodesDomain;iNode++){
          double error = std::abs((myTauBuffer[i][iNode] - tauExpected[i])/tauExpected[i]);
          if(error > 1e-15){
            computeTauWorks = false;
            tauFailed = true;
            if(printTau) {
              std::cout << "In ComputeTauTest, Expected tau[" << i << "]=" << tauExpected[i]
                        << " got tau=" << myTauBuffer[i][iNode] << " at iNode=" << iNode 
                        << " error=" << error << std::endl;
              std::cout << "No further error messages will be generated for this tau index" 
                        << std::endl;
              printTau=false;
            }
          }
        }
      }
      dimComputeStressTensor[numDim-2] = computeTauWorks;
      if(tauFailed){
        std::cout << "Tau failed" << std::endl;
      }
      std::cout << "Done running tests for "<<numDim<<"D ComputeStressTensor" << std::endl;

      std::cout << "Running tests for "<<numDim<<"D ComputeHeatFlux" << std::endl;

      // build a temperature gradient so we can test the calculation of the heat flux vector
      // in 3D...
      // | dT/d(xi) dT/d(eta) dT/d(zeta) |   | [1..1] [2..2] [3..3]
      std::vector<double *> myTempGradBuffer(numDim,NULL);
      std::vector<double> myTempGrad(numDim*numNodesDomain,0);
      for(int i=0;i<numDim;i++){
        myTempGradBuffer[i] = &myTempGrad[i*numNodesDomain];
        double value = static_cast<double>(i+1);
        for(int k=0;k<numNodesDomain;k++){
          myTempGradBuffer[i][k]=value;
        }
      }

      // calculate the heat flux vector
      std::vector<double *> myHeatFluxBuffer(numDim,NULL);
      std::vector<double> myHeatFlux(numDim*numNodesDomain,0);
      for(int i=0;i<numDim;i++){
        myHeatFluxBuffer[i] = &myHeatFlux[i*numNodesDomain];
      }

      if(viscid::util::ComputeHeatFluxBuffer(bufferExtent, numNodes, gridType,gridMetric, 
                                             gridJacobian, Re, myTVBuffers, myTempGradBuffer, 
                                             myHeatFluxBuffer)){
        std::cout << "Calculation of the heat flux vector failed" << std::endl;
      }

      // Construct the expected heat flux vector
      // q_i = mu/Re/Pr * (dT/dx_i) 
      // Re = 100, Pr = 0.75, mu and lambda set above from T=3.5
      // gridmetrics are applied by ComputeHeatFlux, so account for this here
      std::vector<double> heatFluxExpected;
      if (gridType == simulation::grid::RECTILINEAR){
        for(int iDim = 0;iDim < numDim;iDim++){
          double metric = gridMetric[iDim*numNodesDomain];
          double jac    = gridJacobian[0];
          heatFluxExpected.push_back(Cp*tempMu/Pr/Re*(iDim+1)*metric*jac);
        }
      } else if(gridType < simulation::grid::CURVILINEAR) {
        for(int iDim = 0;iDim < numDim;iDim++){
          heatFluxExpected.push_back(Cp*tempMu/Pr/Re*(iDim+1)*gridMetric[iDim]*gridJacobian[0]);
        }
      } else { // must be curvilinear
        for(int iDim = 0;iDim < numDim;iDim++){
          double metric = gridMetric[iDim*numDim*numNodesDomain];
          double jac    = gridJacobian[0];
          heatFluxExpected.push_back(Cp*tempMu/Pr/Re*(iDim+1)*metric*jac);
        }        
      }


      // Compare expected against actual
      bool heatFluxFailed=false;
      bool computeHeatFluxWorks=true;
      for(int i=0;i<numDim;i++){
        bool printHeatFlux=true;
        for(int iNode = 0;iNode < numNodesDomain;iNode++){
          double error = std::abs((myHeatFluxBuffer[i][iNode] - heatFluxExpected[i])/
                                  heatFluxExpected[i]);
          if(error > 1e-15){
            computeHeatFluxWorks = false;
            heatFluxFailed = true;
            if(printHeatFlux) {
              std::cout << "In ComputeHeatFluxTest, Expected heatFlux[" << i << "]=" 
                        << heatFluxExpected[i]
                        << " got heatFlux=" << myHeatFluxBuffer[i][iNode] << " at iNode=" 
                        << iNode << " error=" << error << std::endl;
              std::cout << "No further error messages will be generated for this q index" 
                        << std::endl;
              printHeatFlux=false;
            }
          }
        }
      }
      dimComputeHeatFlux[numDim-2] = computeHeatFluxWorks;
      if(heatFluxFailed){
        std::cout << "Heat flux failed" << std::endl;
      }
      std::cout << "Done running tests for " << numDim 
                << "D ComputeHeatFlux" << std::endl;

      //
      // now test the computation of the viscous flux terms
      // This is completely fabricated and non-physical as the velocity
      // is inconsistent with the velocity gradient.
      //
      std::cout << "Running tests for " << numDim 
                << "D viscous flux" << std::endl;
    
      std::vector<size_t> fortranBufferInterval(2*numDim,1);
      for(int iDim = 0;iDim < numDim;iDim++){
        fortranBufferInterval[2*iDim+1] = bufferInterval[2*iDim+1] + 1;
      }
    
      for(int iDim = 0;iDim < numDim;iDim++){
        vBuffers[iDim] = &velocity[iDim*numNodesDomain];
      }
    
      // New sequence:
      // Compute v_i, tau_i_j and -q_i (done above in original test)
      // Compute the energy flux vector Q_i = (v_j * tau_i_j + q_i)  
      // Use new flux routine (viscidstrongflux1d) to get viscous fluxes:
      // momentum terms = metric * tau
      // energy term = metric * Q
    
      for(int iDim = 0;iDim < numDim;iDim++){
        // Compute q_i += (v_j * tau_i_j)
        for(int iVel = 0;iVel < numDim;iVel++){
          int tensorIndex = iDim*numDim + iVel;
          size_t dirOffset = iVel*numNodesDomain;

          FC_MODULE(operators,ywxpy,OPERATORS,YWXPY)
            (&numDim,&numNodesDomain,&numNodes[0],&fortranBufferInterval[0],
             &velocity[dirOffset],myTauBuffer[tensorIndex],myHeatFluxBuffer[iDim]);
        }
      }

      // Check expected v * tau
      for(int iDim = 0;iDim < numDim;iDim++){
        bool triggered = false;
        double expectedValue = 0;
        for(int iVel = 0;iVel < numDim;iVel++){
          int tensorIndex = iDim*numDim + iVel;
          expectedValue += (vVec[iVel] * tauExpected[tensorIndex]);
        }
        expectedValue += heatFluxExpected[iDim];
        for(size_t iPoint = 0;iPoint < numNodesDomain;iPoint++){
          if((std::abs(myHeatFluxBuffer[iDim][iPoint] - expectedValue)/std::abs(expectedValue)) >
             1e-12){
            dimComputeTEF[nDim-2] = false;
            if(!triggered){
              std::cout << "Failed in " << (iDim+1) << "/" << numDim 
                        << " dimension " << std::endl;
              std::cout << "TEF(" << iPoint << ") [r,e] = [" << myHeatFluxBuffer[iDim][iPoint]
                        << "," << expectedValue << "]" << std::endl;
              triggered = true;
            }
          }
        }
      }
    
    
      for(int iDim = 0;iDim < numDim;iDim++){

        double metricX = gridMetric[iDim];

        if(gridType == simulation::grid::RECTILINEAR) {
          metricX = gridMetric[iDim*numNodesDomain];
        } else if (gridType == simulation::grid::CURVILINEAR) {
          metricX = gridMetric[iDim*numDim*numNodesDomain];
        }

        std::vector<double> viscidFlux(numNodesDomain*(numDim+2),0.0);
        std::vector<double> expectedFlux((numDim+2)*numNodesDomain,0.0);
        std::vector<double> scaledTau(numNodesDomain,0.0);

        for(size_t iPoint = 0;iPoint < numNodesDomain;iPoint++){
          size_t fOffset = 0;

          for(int ivDim = 0;ivDim < numDim;ivDim++){
            size_t iOffset = ivDim*numNodesDomain;
            fOffset = iOffset + numNodesDomain;
            size_t vectorOffset = ivDim*numNodesDomain;

            // momentum
            int tauIndex=iDim*numDim+ivDim;
            expectedFlux[iPoint+fOffset] += tauExpected[tauIndex]*metricX;

            // energy
            expectedFlux[iPoint+(numDim+1)*numNodesDomain] += 
              tauExpected[tauIndex]*velocity[vectorOffset+iPoint]*metricX;
          }

          // thermal conduction
          expectedFlux[iPoint+(numDim+1)*numNodesDomain] +=
            heatFluxExpected[iDim]*metricX;

          fOffset += numNodesDomain;
        }

        int fluxDim = iDim + 1;
        int tau1=iDim*numDim;
        int tau2=tau1+1;
        int tau3=tau2+1;
      
        FC_MODULE(viscid,strongflux1d,VISCID,STRONGFLUX1D)
          (&numDim,&fluxDim,&numNodes[0],&numNodesDomain,&fortranBufferInterval[0],
           &gridType,&gridMetric[0],&myTau[0],&myHeatFlux[0],&viscidFlux[0]);

        std::vector<bool> uniformFlux(numDim+2,true);
        for(int fDim = 0 ;fDim < (numDim+2);fDim++){
          size_t fOffset = fDim * numNodesDomain;
          bool triggered = false;
          for(size_t iPoint = 0;iPoint < numNodesDomain;iPoint++){
            double error= std::abs((viscidFlux[fOffset+iPoint] - expectedFlux[fOffset+iPoint])/
                                   expectedFlux[fOffset+iPoint]);
            if(error > 1e-10){
              if(!triggered ){
                std::cout << "iDim " << iDim << " Flux Dimension " << fDim << std::endl;
                std::cout << std::scientific << std::setprecision(16)
                          << "\tActual:     " << viscidFlux[fOffset+iPoint] << std::endl
                          << "\tExpected:   " << expectedFlux[fOffset+iPoint] << std::endl
                          << "\tError:      " << error << std::endl; 
                triggered = true;
              }
              uniformFlux[fDim] = false;
            }
          }
        }
        for(int fDim = 0;fDim < (numDim+2);fDim++){
          if(!uniformFlux[fDim]){
            dimViscidUniformFlux[numDim-2] = false;
            std::cout << "uniform flux failed in dim,fdim " << iDim 
                      << "," << fDim << std::endl;
          }
        }
        std::cout << "Done running tests for "<<numDim<<"D viscous flux" << std::endl;

      } // loop over iDim

    } // nDim

    bool computeDVWorks = true;
    bool computeTVWorks = true;
    bool computeStressTensorWorks = true;
    bool computeHeatFluxWorks = true;
    bool viscidUniformFluxWorks = true;
    bool computeTotalEnergyFlux = true;
    // Make sure these tests worked in 2d and 3d
    for(int nDim = 2;nDim <= 3;nDim++){
      // compute DV tests are temporay until I have computeTV working
      computeDVWorks = computeDVWorks && dimComputeDV2[nDim-2];
      computeTVWorks = computeTVWorks && dimComputeTV[nDim-2];
      computeTotalEnergyFlux = computeTotalEnergyFlux && dimComputeTEF[nDim-2];
      computeStressTensorWorks = computeStressTensorWorks && dimComputeStressTensor[nDim-2];
      computeHeatFluxWorks = computeHeatFluxWorks && dimComputeHeatFlux[nDim-2];
      viscidUniformFluxWorks = viscidUniformFluxWorks && dimViscidUniformFlux[nDim-2];
    }
    if(computeStressTensorWorks && computeHeatFluxWorks && viscidUniformFluxWorks)
      metricSuccess[gridType] = true;
    else {
      metricSuccess[gridType] = false;
      testResult = false;
    }

  } // gridType
  
  
  serialUnitResults.UpdateResult("Viscid:Kernels:GridMetrics",testResult);

//   serialUnitResults.UpdateResult("Viscid:Kernels:HeatFlux:GridMetrics",computeHeatFluxWorks);
//   serialUnitResults.UpdateResult("Viscid:Kernels:UniformViscidFlux:GridMetrics",viscidUniformFluxWorks);

};

void TestVelocityGradient(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{
  std::cout << "TestVelocityGradient..." << std::endl;

  typedef plascom2::operators::sbp::base  operator_t;
  typedef simulation::grid::halo          halo_t;
  typedef simulation::state::base         state_t;
  typedef pcpp::ParallelGlobalType        global_t;
  typedef pcpp::CommunicatorType          comm_t;
  typedef pcpp::IndexIntervalType         interval_t;

  typedef simulation::grid::parallel_blockstructured         grid_t;
  typedef simulation::domain::base<grid_t,state_t>  domain_t;
  typedef euler::rhs<grid_t,state_t,operator_t>              rhs_t;

  std::ostringstream messageStream;

  grid_t     simGrid2d;
  grid_t     simGrid3d;
  operator_t sbpOperator2d;
  operator_t sbpOperator3d;
  halo_t     &simHalo2d(simGrid2d.Halo());
  halo_t     &simHalo3d(simGrid3d.Halo());

  state_t    simState2d;
  state_t    simState3d;
  state_t    paramState2d;
  state_t    paramState3d;
  
  bool parTestViscidRHS = true;
  bool VelocityGradientWorks = true;
  bool TemperatureGradientWorks = true;

  std::cout << "Starting 2D linear Velocity Gradient Test" << std::endl;

  // sets up a 2D grid with 4 points in each direction
  int numNodes = 21;
  std::vector<size_t> gridSizes(2,numNodes);
  simGrid2d.SetGridSizes(gridSizes);

  std::vector<double> physicalExtent(4,1.0);
  physicalExtent[0] = 0.0;
  physicalExtent[2] = 0.0;
  simGrid2d.SetPhysicalExtent(physicalExtent);

  std::vector<bool> periodicDirs(2,false);
  simGrid2d.SetPeriodicDirs(periodicDirs);

  if(euler::util::InitializeSimulationFixtures(simGrid2d,sbpOperator2d,2,testComm,&messageStream)){
    parTestViscidRHS = false;
    std::cout << "Failed to setup 2D grid in TestViscidRHS" << std::endl;
    return;
  }

  std::cout << messageStream.str() << std::endl;
  messageStream.str("");
  messageStream.clear();

  const pcpp::IndexIntervalType &partitionBufferInterval2d(simGrid2d.PartitionBufferInterval());
  const pcpp::IndexIntervalType &partitionInterval2d(simGrid2d.PartitionInterval());
  int numDim = gridSizes.size();


  if(testfixtures::viscid::SetupViscidState(simGrid2d,simState2d,paramState2d,0)){
    parTestViscidRHS = false;
    std::cout << "Failed to setup state in TestViscidRHS" << std::endl;
    return;
  }
  
  size_t numPointsBuffer = simGrid2d.BufferSize();
  const std::vector<size_t> &bufferSizes(simGrid2d.BufferSizes());

  double t = 0.0;
  double gamma = 1.4;
  double dt = .001;
  double cfl = 1.0;
  double Re = 100;
  double power = 0.666;
  double beta = 5;
  double bulkViscFac = 2;
  double Pr = 0.75;
  std::vector<double> &dX(simGrid2d.DX());

  // Add some fields that are not usual
  paramState2d.AddField("dX",'s',numDim,8,"space");
  paramState2d.Create(numPointsBuffer,0);

  paramState2d.SetFieldBuffer("gamma",&gamma);
  paramState2d.SetFieldBuffer("inputDT",&dt);
  paramState2d.SetFieldBuffer("inputCFL",&cfl);
  paramState2d.SetFieldBuffer("dX",&dX[0]);
  paramState2d.SetFieldBuffer("refRe",&Re);
  paramState2d.SetFieldBuffer("refPr",&Pr);
  paramState2d.SetFieldBuffer("power",&power);
  paramState2d.SetFieldBuffer("beta",&beta);
  paramState2d.SetFieldBuffer("bulkViscFac",&bulkViscFac);

  std::vector<double> rho(numPointsBuffer,1.0);
  std::vector<double> rhoV(numDim*numPointsBuffer,0.0);
  std::vector<double> rhoE(numPointsBuffer,1.0/(gamma-1.0));

  std::vector<double> T(numPointsBuffer,0.0);
  std::vector<double> p(numPointsBuffer,0.0);
  std::vector<double> V(3*numPointsBuffer,0.0);
  std::vector<double> rhom1(numPointsBuffer,0.0);

  simState2d.SetFieldBuffer("simTime",&t);
  simState2d.SetFieldBuffer("rho",&rho[0]);
  simState2d.SetFieldBuffer("rhoV",&rhoV[0]);
  simState2d.SetFieldBuffer("rhoE",&rhoE[0]);

  simState2d.SetFieldBuffer("pressure",&p[0]);
  simState2d.SetFieldBuffer("temperature",&T[0]);
  simState2d.SetFieldBuffer("velocity",&V[0]);
  simState2d.SetFieldBuffer("rhom1",&rhom1[0]);

  state_t rhsState(simState2d);


  // -- Set up the RHS --
  rhs_t viscidRHS2d;
  
  if(viscidRHS2d.Initialize(simGrid2d,simState2d,paramState2d,sbpOperator2d)){
    parTestViscidRHS = false;
    std::cout << "Failed to initialize RHS in TestViscidRHS" << std::endl;
    return;
  }

  // Insert thread partitioning here 
  // - initialize RHS and grid thread intervals
  std::vector<int> numThreadsDim;
  simGrid2d.SetupThreads(numThreadsDim);
  viscidRHS2d.InitThreadIntervals();
 
  if(viscidRHS2d.SetupRHS(simGrid2d,simState2d,rhsState)){
    parTestViscidRHS = false;
    std::cout << "Failed to setup RHS in TestViscidRHS" << std::endl;
    return;
  }

  parallelUnitResults.UpdateResult("Viscid:RHS:RunsInParallel",parTestViscidRHS);

  // 
  // setup a 2D linear velocity field so we can test the velocity gradient
  // setup a 2D linear temperature field so we can test the temperature gradient
  // leave the density alone, update the energy to be consistent with the temperature
  //
  std::vector<double> velGradExpected(numDim*numDim*numPointsBuffer,-999.0);
  std::vector<double> tempGradExpected(numDim*numPointsBuffer,-999.0);
  std::vector<double> uAmp; 
  std::vector<double> vAmp;
  std::vector<double> tempAmp;
  uAmp.push_back(1.0);
  uAmp.push_back(2.0);
  uAmp.push_back(0.0);
  vAmp.push_back(3.0);
  vAmp.push_back(4.0);
  vAmp.push_back(0.0);
  tempAmp.push_back(5.0);
  tempAmp.push_back(6.0);
  tempAmp.push_back(0.0);
  double x0 = 0.;
  double y0 = 0.;
  size_t iStart = partitionBufferInterval2d[0].first;
  size_t iEnd   = partitionBufferInterval2d[0].second;
  size_t jStart = partitionBufferInterval2d[1].first;
  size_t jEnd   = partitionBufferInterval2d[1].second;
  for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
    size_t jBufferIndex = jIndex*bufferSizes[0];
    size_t jPartIndex = (jIndex - jStart) + partitionInterval2d[1].first;
    double y = jPartIndex*dX[1];
    for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
      size_t bufferIndex = jBufferIndex + iIndex;
      size_t iPartIndex = (iIndex - iStart) + partitionInterval2d[0].first;
      double x = iPartIndex*dX[0];

      rhoV[bufferIndex] = rho[bufferIndex]*
        testfixtures::viscid::Linear(uAmp,x0,y0,0.,x,y,0.,-1);
      rhoV[bufferIndex+numPointsBuffer] = rho[bufferIndex]*
        testfixtures::viscid::Linear(vAmp,x0,y0,0.,x,y,0.,-1);
      // just the internal part, add a linear offset
      rhoE[bufferIndex] += rho[bufferIndex]/(gamma-1)*
        testfixtures::viscid::Linear(tempAmp,x0,y0,0.,x,y,0.,-1);

      //std::cout << "Expected velocity (" << rhoV[bufferIndex] << "," 
                //<< rhoV[bufferIndex+numPointsBuffer] << "," 
                //<< rhoV[bufferIndex+2*numPointsBuffer] 
                //<< ") at (" << x << "," << y << "," << z << ")" << std::endl;
      rhoE[bufferIndex] = rhoE[bufferIndex] + (rhoV[bufferIndex]*rhoV[bufferIndex] +
                             rhoV[bufferIndex+numPointsBuffer]*
                             rhoV[bufferIndex+numPointsBuffer])/(2.0*rho[bufferIndex]);

      // calculate the exact velocity gradient for comparison
      // one big flat vector with numPointsBuffer*numDin^2 size
      // each {} is numPoints in size
      // [ {du/dx} {dv/dx} {du/dy} {dv/dy} ]
      // we remove the grid space as the metrics are applied inside ComputeStressTensor
      for(int dIndex = 0; dIndex<numDim; dIndex++){
        velGradExpected[bufferIndex+dIndex*numDim*numPointsBuffer] = 
          testfixtures::viscid::Linear(uAmp,x0,y0,0.,x,y,0.,dIndex)*dX[dIndex];
        velGradExpected[bufferIndex+numPointsBuffer+dIndex*numDim*numPointsBuffer] = 
          testfixtures::viscid::Linear(vAmp,x0,y0,0.,x,y,0.,dIndex)*dX[dIndex];
      }

      // calculate the exact temperature gradient for comparison
      // one big flat vector with numPointsBuffer*numDim size
      // each {} is numPoints in size
      // [ {dT/dx} {dT/dy} ]
      // we remove the grid space as the metrics are applied inside ComputeStressTensor
      for(int dIndex = 0; dIndex<numDim; dIndex++){
        tempGradExpected[bufferIndex+dIndex*numPointsBuffer] = 
          testfixtures::viscid::Linear(tempAmp,x0,y0,0.,x,y,0.,dIndex)*dX[dIndex];
      }
    }
  }

  int returnCode = viscidRHS2d.ComputeDV(0);
  if(returnCode > 0){
    VelocityGradientWorks = false;
    std::cout << "ComputeVelGrad failed when computing DV " << std::endl;
    return;
  }

  std::cout << "Calculating Velocity Gradient" << std::endl;
  returnCode = viscidRHS2d.ComputeVelGrad(0);
  if(returnCode > 0){
    VelocityGradientWorks = false;
    return;
  }

  VelocityGradientWorks = CheckResult(VelocityGradientWorks,testComm);
  if(!VelocityGradientWorks){
    std::cout << "Velocity Gradient Calculation Failed" << std::endl;
  }


  // check the velocity gradient
  // first get the gradient from the rhs object
  // then check against the expected answer
  std::cout << "Checking Velocity Gradient calculation" << std::endl;
  const std::vector<double *> &velGradBuffers(viscidRHS2d.VelGradBuffers());

  if(VelocityGradientWorks){
    bool printVelGrad = true;
    for(int iDim = 0;iDim < numDim;iDim++){
      const double* velGradData = velGradBuffers[iDim];
      
      for(int jDim = 0;jDim < numDim;jDim++){
        const double* velGradDataComponent = &velGradData[jDim*numPointsBuffer];
        
        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jBufferIndex = jIndex*bufferSizes[0];
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jBufferIndex + iIndex;
            size_t expectedIndex = numPointsBuffer*iDim*numDim+jDim*numPointsBuffer+bufferIndex;
            
            double error = std::abs(velGradDataComponent[bufferIndex] - 
                                    velGradExpected[expectedIndex])/velGradExpected[expectedIndex];
            
            if(error > 1e-12){
              VelocityGradientWorks = false;
              if(printVelGrad) {
                std::cout << std::scientific << std::setprecision(16)
                          << "In ComputeVelGrad: "<< std::endl
                          << "Expected velGrad=" << velGradExpected[expectedIndex] << std::endl
                          << " got velGrad=" << velGradDataComponent[bufferIndex] << std::endl
                          << " at (i,j,bufferIndex)=(" << iDim << "," << jDim << "," 
                          << bufferIndex << ")" << std::endl;
                std::cout << "Error =" << error << std::endl;;
                std::cout << "No further error messages will be generated for velGrad 2D Linear" << std::endl;
                printVelGrad=false;
              }
            }
          }
        }
      }
    }
  }
  VelocityGradientWorks = CheckResult(VelocityGradientWorks,testComm);
  if(!VelocityGradientWorks){
    std::cout << "Velocity Gradient Compare Failed" << std::endl;
  }

  // Compute the temperature gradient
  std::cout << "Calculating Temperature Gradient" << std::endl;
  returnCode = viscidRHS2d.ComputeTemperatureGrad(0);
  if(returnCode > 0){
    TemperatureGradientWorks = false;
    return;
  }

  TemperatureGradientWorks = CheckResult(TemperatureGradientWorks,testComm);
  if(!TemperatureGradientWorks){
    std::cout << "Temperature Gradient Calculation Failed" << std::endl;
  }

  // check the temperature gradient
  // first get the gradient from the rhs object
  // then check against the expected answer
  std::cout << "Checking Temperature Gradient calculation" << std::endl;
  const std::vector<double *> &tempGradBuffers(viscidRHS2d.TemperatureGradBuffers());

  if(TemperatureGradientWorks){
    bool printTempGrad = true;
    for(int iDim = 0;iDim < numDim;iDim++){
      const double* tempGradDataComponent = tempGradBuffers[iDim];
      
      for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
        size_t jBufferIndex = jIndex*bufferSizes[0];
        for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
          size_t bufferIndex = jBufferIndex + iIndex;
          size_t expectedIndex = bufferIndex+numPointsBuffer*iDim;
          
          double error = std::abs(tempGradDataComponent[bufferIndex] - 
            tempGradExpected[expectedIndex])/tempGradExpected[expectedIndex];
          
          if(error > 1e-12){
            TemperatureGradientWorks = false;
            if(printTempGrad) {
              std::cout << std::scientific << std::setprecision(16)
                        << "In ComputeTemperatureGrad: "<< std::endl
                        << "Expected tempGrad=" << tempGradExpected[expectedIndex] << std::endl
                        << " got tempGrad=" << tempGradDataComponent[bufferIndex] << std::endl
                        << " at (i,bufferIndex)=(" << iDim << "," 
                        << bufferIndex << ")" << std::endl;
              std::cout << "Temperature " << T[bufferIndex] << std::endl;
              std::cout << "Error =" << error << std::endl;;
              std::cout << "No further error messages will be generated for tempGrad 2D Linear" 
                        << std::endl;
              printTempGrad=false;
            }
          }
        }
      }
    }
  }
  TemperatureGradientWorks = CheckResult(TemperatureGradientWorks,testComm);
  if(!TemperatureGradientWorks){
    std::cout << "Temperature Gradient Compare Failed" << std::endl;
  }
  
  std::cout << "Done with 2D linear Velocity Gradient Test" << std::endl;
  std::cout << "Starting 3D linear Velocity Gradient Test" << std::endl;

  numNodes = 21;
  gridSizes.resize(3,numNodes);
  simGrid3d.SetGridSizes(gridSizes);

  physicalExtent.resize(6,1.0);
  physicalExtent[0] = 0.0;
  physicalExtent[2] = 0.0;
  physicalExtent[4] = 0.0;
  simGrid3d.SetPhysicalExtent(physicalExtent);

  periodicDirs.resize(3,false);
  simGrid3d.SetPeriodicDirs(periodicDirs);

  if(euler::util::InitializeSimulationFixtures(simGrid3d,sbpOperator3d,2,testComm,&messageStream)){
    parTestViscidRHS = false;
    std::cout << "Failed to setup 3D grid in TestViscidRHS" << std::endl;
    return;
  }

  std::cout << messageStream.str() << std::endl;
  messageStream.str("");
  messageStream.clear();

  const pcpp::IndexIntervalType &partitionBufferInterval3d(simGrid3d.PartitionBufferInterval());
  const pcpp::IndexIntervalType &partitionInterval3d(simGrid3d.PartitionInterval());
  numDim = gridSizes.size();

  if(testfixtures::viscid::SetupViscidState(simGrid3d,simState3d,paramState3d,0)){
    parTestViscidRHS = false;
    std::cout << "Failed to setup state in TestViscidRHS" << std::endl;
    return;
  }

  
  numPointsBuffer = simGrid3d.BufferSize();
  //bufferSizes.resize(simGrid3d.BufferSizes());
  const std::vector<size_t> &bufferSizes3d(simGrid3d.BufferSizes());

  //std::vector<double> &dX.resize(simGrid3d.DX());
  std::vector<double> &dX3d(simGrid3d.DX());

  // Add some fields that are not usual
  paramState3d.AddField("dX",'s',3,8,"space");
  paramState3d.Create(numPointsBuffer,0);

  paramState3d.SetFieldBuffer("gamma",&gamma);
  paramState3d.SetFieldBuffer("inputDT",&dt);
  paramState3d.SetFieldBuffer("inputCFL",&cfl);
  paramState3d.SetFieldBuffer("dX",&dX3d[0]);
  paramState3d.SetFieldBuffer("refRe",&Re);
  paramState3d.SetFieldBuffer("refPr",&Pr);
  paramState3d.SetFieldBuffer("power",&power);
  paramState3d.SetFieldBuffer("beta",&beta);
  paramState3d.SetFieldBuffer("bulkViscFac",&bulkViscFac);

  rho.resize(numPointsBuffer,1.0);
  rhoV.resize(3*numPointsBuffer,0.0);
  rhoE.resize(numPointsBuffer);
  rhoE.assign(rhoE.size(),1.0/(gamma-1));

  T.resize(numPointsBuffer,0.0);
  p.resize(numPointsBuffer,0.0);
  V.resize(3*numPointsBuffer,0.0);

  rhom1.resize(numPointsBuffer,0.0);

  simState3d.SetFieldBuffer("simTime",&t);
  simState3d.SetFieldBuffer("rho",&rho[0]);
  simState3d.SetFieldBuffer("rhoV",&rhoV[0]);
  simState3d.SetFieldBuffer("rhoE",&rhoE[0]);

  simState3d.SetFieldBuffer("pressure",&p[0]);
  simState3d.SetFieldBuffer("temperature",&T[0]);
  simState3d.SetFieldBuffer("velocity",&V[0]);
  simState3d.SetFieldBuffer("rhom1",&rhom1[0]);

  state_t rhsState3d(simState3d);


  // -- Set up the RHS --
  rhs_t viscidRHS3d;
  
  if(viscidRHS3d.Initialize(simGrid3d,simState3d,paramState3d,sbpOperator3d)){
    parTestViscidRHS = false;
    std::cout << "Failed to initialize RHS in TestViscidRHS" << std::endl;
    return;
  }

  std::vector<int> numThreadsDim3;
  simGrid3d.SetupThreads(numThreadsDim3);
  viscidRHS3d.InitThreadIntervals();

  if(viscidRHS3d.SetupRHS(simGrid3d,simState3d,rhsState3d)){
    parTestViscidRHS = false;
    std::cout << "Failed to setup RHS in TestViscidRHS" << std::endl;
    return;
  }

  parallelUnitResults.UpdateResult("Viscid:RHS:RunsInParallel",parTestViscidRHS);

  // 
  // setup a 3D linear velocity field so we can test the velocity gradient
  // setup a 3D linear temperature field so we can test the temperature gradient
  // leave the density alone, update the energy
  //
  velGradExpected.resize(numDim*numDim*numPointsBuffer,-999.0);
  tempGradExpected.resize(numDim*numPointsBuffer,-999.0);
  uAmp.resize(0); 
  vAmp.resize(0);
  tempAmp.resize(0);
  std::vector<double> wAmp;
  uAmp.push_back(1.0);
  uAmp.push_back(2.0);
  uAmp.push_back(3.0);
  vAmp.push_back(4.0);
  vAmp.push_back(5.0);
  vAmp.push_back(6.0);
  wAmp.push_back(7.0);
  wAmp.push_back(8.0);
  wAmp.push_back(9.0);
  tempAmp.push_back(10.0);
  tempAmp.push_back(11.0);
  tempAmp.push_back(12.0);
  x0 = 0.;
  y0 = 0.;
  double z0 = 0.;
  iStart = partitionBufferInterval3d[0].first;
  iEnd   = partitionBufferInterval3d[0].second;
  jStart = partitionBufferInterval3d[1].first;
  jEnd   = partitionBufferInterval3d[1].second;
  size_t kStart = partitionBufferInterval3d[2].first;
  size_t kEnd   = partitionBufferInterval3d[2].second;
  size_t nPlane = bufferSizes3d[0]*bufferSizes3d[1];
  for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
    size_t kBufferIndex = kIndex*nPlane;
    size_t kPartIndex = (kIndex - kStart) + partitionInterval3d[2].first;
    double z = kPartIndex*dX3d[2]; 
    for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
      size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes3d[0];
      size_t jPartIndex = (jIndex - jStart) + partitionInterval3d[1].first;
      double y = jPartIndex*dX3d[1];
      for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
        size_t bufferIndex = jkBufferIndex + iIndex;
        size_t iPartIndex = (iIndex - iStart) + partitionInterval3d[0].first;
        double x = iPartIndex*dX3d[0];

        rhoV[bufferIndex] = rho[bufferIndex]*
          testfixtures::viscid::Linear(uAmp,x0,y0,z0,x,y,z,-1);
        rhoV[bufferIndex+numPointsBuffer] = rho[bufferIndex]*
          testfixtures::viscid::Linear(vAmp,x0,y0,z0,x,y,z,-1);
        rhoV[bufferIndex+2*numPointsBuffer] = rho[bufferIndex]*
          testfixtures::viscid::Linear(wAmp,x0,y0,z0,x,y,z,-1);

        // just the internal part, add a linear offset
        rhoE[bufferIndex] += rho[bufferIndex]/(gamma-1)*
          testfixtures::viscid::Linear(tempAmp,x0,y0,z0,x,y,z,-1);

        rhoE[bufferIndex] = rhoE[bufferIndex] + (rhoV[bufferIndex]*rhoV[bufferIndex] +
                               rhoV[bufferIndex+numPointsBuffer]*
                               rhoV[bufferIndex+numPointsBuffer] +
                               rhoV[bufferIndex+2*numPointsBuffer]*
                               rhoV[bufferIndex+2*numPointsBuffer])/(2.0*rho[bufferIndex]);

        // calculate the exact velocity gradient for comparison
        // one big flat vector with numPointsBuffer*numDin^2 size
        // each {} is numPoints in size
        // [ {du/dx} {dv/dx} {dw/dx} {du/dy} {dv/dy} {dw/dy} {du/dz} {dv/dz} {dw/dz} ]
        for(int dIndex = 0; dIndex<numDim; dIndex++){
          velGradExpected[bufferIndex+dIndex*numDim*numPointsBuffer] = 
            testfixtures::viscid::Linear(uAmp,x0,y0,z0,x,y,z,dIndex)*dX3d[dIndex];
          velGradExpected[bufferIndex+numPointsBuffer+dIndex*numDim*numPointsBuffer] = 
            testfixtures::viscid::Linear(vAmp,x0,y0,z0,x,y,z,dIndex)*dX3d[dIndex];
          velGradExpected[bufferIndex+2*numPointsBuffer+dIndex*numDim*numPointsBuffer] = 
            testfixtures::viscid::Linear(wAmp,x0,y0,z0,x,y,z,dIndex)*dX3d[dIndex];
        }

        // calculate the exact temperature gradient for comparison
        // one big flat vector with numPointsBuffer*numDim size
        // each {} is numPoints in size
        // [ {dT/dx} {dT/dy} {dT/dz}]
        // we remove the grid space as the metrics are applied inside ComputeStressTensor
        for(int dIndex = 0; dIndex<numDim; dIndex++){
          tempGradExpected[bufferIndex+dIndex*numPointsBuffer] = 
            testfixtures::viscid::Linear(tempAmp,x0,y0,z0,x,y,z,dIndex)*dX3d[dIndex];
        }
      }
    }
  }

  returnCode = viscidRHS3d.ComputeDV(0);
  if(returnCode > 0){
    VelocityGradientWorks = false;
    std::cout << "ComputeVelGrad failed when computing DV " << std::endl;
    return;
  }

  std::cout << "Calculating Velocity Gradient" << std::endl;
  returnCode = viscidRHS3d.ComputeVelGrad(0);
  if(returnCode > 0){
    VelocityGradientWorks = false;
    return;
  }

  VelocityGradientWorks = CheckResult(VelocityGradientWorks,testComm);
  if(!VelocityGradientWorks){
    std::cout << "Velocity Gradient Calculation Failed" << std::endl;
  }

  if(VelocityGradientWorks){
    std::cout << "Checking Velocity Gradient calculation" << std::endl;
    // check the velocity gradient
    // first get the gradient from the rhs object
    // then check against the expected answer
    const std::vector<double *> &velGradBuffers3d(viscidRHS3d.VelGradBuffers());
    
    bool printVelGrad = true;
    for(int iDim = 0;iDim < numDim;iDim++){
      const double* velGradData3d = velGradBuffers3d[iDim];
      
      for(int jDim = 0;jDim < numDim;jDim++){
        const double* velGradDataComponent3d = &velGradData3d[jDim*numPointsBuffer];
        
        for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
          size_t kBufferIndex = kIndex*nPlane;
          for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
            size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes3d[0];
            for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
              size_t bufferIndex = jkBufferIndex + iIndex;
              size_t expectedIndex = numPointsBuffer*iDim*numDim+jDim*numPointsBuffer+bufferIndex;
              
              double error = std::abs(velGradDataComponent3d[bufferIndex] - 
                velGradExpected[expectedIndex])/velGradExpected[expectedIndex];
              
              if(error > 1e-12){
                VelocityGradientWorks = false;
                if(printVelGrad) {
                  std::cout << std::scientific << std::setprecision(16)
                            << "In ComputeVelGrad: "<< std::endl
                            << "Expected velGrad=" << velGradExpected[expectedIndex] << std::endl
                            << " got velGrad=" << velGradDataComponent3d[bufferIndex] << std::endl
                            << " at (i,j,bufferIndex)=(" << iDim << "," << jDim << "," 
                            << bufferIndex << ")" << std::endl;
                  std::cout << "Error =" << error << std::endl;;
                  std::cout << "No further error messages will be generated for velGrad 3d" << std::endl;
                  printVelGrad=false;
                }
              }
            }
          }
        }
      }
    }
  }

  // Compute the temperature gradient
  std::cout << "Calculating Temperature Gradient" << std::endl;
  returnCode = viscidRHS3d.ComputeTemperatureGrad(0);
  if(returnCode > 0){
    TemperatureGradientWorks = false;
    return;
  }

  TemperatureGradientWorks = CheckResult(TemperatureGradientWorks,testComm);
  if(!TemperatureGradientWorks){
    std::cout << "Temperature Gradient Calculation Failed" << std::endl;
  }

  // check the temperature gradient
  // first get the gradient from the rhs object
  // then check against the expected answer
  std::cout << "Checking Temperature Gradient calculation" << std::endl;
  const std::vector<double *> &tempGradBuffers3d(viscidRHS3d.TemperatureGradBuffers());

  if(TemperatureGradientWorks){
    bool printTempGrad = true;
    for(int iDim = 0;iDim < numDim;iDim++){
      const double* tempGradDataComponent = tempGradBuffers3d[iDim];

      for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
        size_t kBufferIndex = kIndex*nPlane;
        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes3d[0];
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jkBufferIndex + iIndex;
            size_t expectedIndex = bufferIndex + numPointsBuffer*iDim;
    
            double error = std::abs(tempGradDataComponent[bufferIndex] - 
              tempGradExpected[expectedIndex])/tempGradExpected[expectedIndex];
          
            if(error > 1e-12){
              TemperatureGradientWorks = false;
              if(printTempGrad) {
                std::cout << std::scientific << std::setprecision(16)
                          << "In ComputeTemperatureGrad: "<< std::endl
                          << "Expected tempGrad=" << tempGradExpected[expectedIndex] << std::endl
                          << " got tempGrad=" << tempGradDataComponent[bufferIndex] << std::endl
                          << " at (i,bufferIndex)=(" << iDim << "," 
                          << bufferIndex << ")" << std::endl;
                std::cout << "Temperature " << T[bufferIndex] << std::endl;
                std::cout << "Error =" << error << std::endl;;
                std::cout << "No further error messages will be generated for tempGrad 3D Linear" 
                          << std::endl;
                printTempGrad=false;
              }
            }
          }
        }
      }
    }
  }

  std::cout << "Done with 3D linear Velocity Gradient Test" << std::endl;

  TemperatureGradientWorks = CheckResult(TemperatureGradientWorks,testComm);
  if(!TemperatureGradientWorks){
    std::cout << "Temperature Gradient Compare Failed" << std::endl;
  }

  VelocityGradientWorks = CheckResult(VelocityGradientWorks,testComm);
  if(!VelocityGradientWorks){
    std::cout << "Velocity Gradient Compare Failed" << std::endl;
  }
  

  parallelUnitResults.UpdateResult("Viscid:RHS:VelocityGradient",VelocityGradientWorks);
  parallelUnitResults.UpdateResult("Viscid:RHS:TemperatureGradient",TemperatureGradientWorks);

};

void TestVelocityGradientPeriodic(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{
  std::cout << "TestVelocityGradientPeriodic..." << std::endl;

  typedef plascom2::operators::sbp::base  operator_t;
  typedef simulation::grid::halo          halo_t;
  typedef simulation::state::base         state_t;
  typedef pcpp::ParallelGlobalType        global_t;
  typedef pcpp::CommunicatorType          comm_t;
  typedef pcpp::IndexIntervalType         interval_t;

  typedef simulation::grid::parallel_blockstructured         grid_t;
  typedef simulation::domain::base<grid_t,state_t>  domain_t;
  typedef euler::rhs<grid_t,state_t,operator_t>              rhs_t;

  std::ostringstream messageStream;

  grid_t     simGrid2d;
  grid_t     simGrid3d;
  operator_t sbpOperator2d;
  operator_t sbpOperator3d;
  halo_t     &simHalo2d(simGrid2d.Halo());
  halo_t     &simHalo3d(simGrid3d.Halo());

  state_t    simState2d;
  state_t    simState3d;
  state_t    paramState2d;
  state_t    paramState3d;
  
  bool parTestViscidRHS = true;
  bool VelocityGradientWorks = true;
  bool TemperatureGradientWorks = true;

  std::cout << "Starting 2D periodic Velocity Gradient Test" << std::endl;

  // sets up a 2D grid with 4 points in each direction
  int numNodes = 31;
  std::vector<size_t> gridSizes(2,numNodes);
  simGrid2d.SetGridSizes(gridSizes);

  std::vector<double> physicalExtent(4,1.0);
  physicalExtent[0] = 0.0;
  physicalExtent[2] = 0.0;
  simGrid2d.SetPhysicalExtent(physicalExtent);

  std::vector<bool> periodicDirs(2,true);
  simGrid2d.SetPeriodicDirs(periodicDirs);
   
  if(euler::util::InitializeSimulationFixtures(simGrid2d,sbpOperator2d,2,testComm,&messageStream)){
    parTestViscidRHS = false;
    std::cout << "Failed to setup 2D grid in TestViscidRHS" << std::endl;
    return;
  }

  std::cout << messageStream.str() << std::endl;
  messageStream.str("");
  messageStream.clear();

  const pcpp::IndexIntervalType &partitionBufferInterval2d(simGrid2d.PartitionBufferInterval());
  const pcpp::IndexIntervalType &partitionInterval2d(simGrid2d.PartitionInterval());
  int numDim = gridSizes.size();

  if(testfixtures::viscid::SetupViscidState(simGrid2d,simState2d,paramState2d,0)){
    parTestViscidRHS = false;
    std::cout << "Failed to setup state in TestViscidRHS" << std::endl;
    return;
  }
  
  size_t numPointsBuffer = simGrid2d.BufferSize();
  const std::vector<size_t> &bufferSizes(simGrid2d.BufferSizes());

  double t = 0.0;
  double gamma = 1.4;
  double dt = .001;
  double cfl = 1.0;
  double Re = 100;
  double Pr = 0.75;
  double power = 0.666;
  double beta = 5;
  double bulkViscFac = 2;
  std::vector<double> &dX(simGrid2d.DX());

  // Add some fields that are not usual
  paramState2d.AddField("dX",'s',numDim,8,"space");
  paramState2d.Create(numPointsBuffer,0);

  paramState2d.SetFieldBuffer("gamma",&gamma);
  paramState2d.SetFieldBuffer("inputDT",&dt);
  paramState2d.SetFieldBuffer("inputCFL",&cfl);
  paramState2d.SetFieldBuffer("dX",&dX[0]);
  paramState2d.SetFieldBuffer("refRe",&Re);
  paramState2d.SetFieldBuffer("refPr",&Pr);
  paramState2d.SetFieldBuffer("power",&power);
  paramState2d.SetFieldBuffer("beta",&beta);
  paramState2d.SetFieldBuffer("bulkViscFac",&bulkViscFac);

  std::vector<double> rho(numPointsBuffer,1.0);
  std::vector<double> rhoV(numDim*numPointsBuffer,0.0);
  std::vector<double> rhoE(numPointsBuffer,1.0/(gamma-1.0));

  std::vector<double> T(numPointsBuffer,0.0);
  std::vector<double> p(numPointsBuffer,0.0);
  std::vector<double> V(numDim*numPointsBuffer,0.0);
  std::vector<double> rhom1(numPointsBuffer,0.0);

  simState2d.SetFieldBuffer("simTime",&t);
  simState2d.SetFieldBuffer("rho",&rho[0]);
  simState2d.SetFieldBuffer("rhoV",&rhoV[0]);
  simState2d.SetFieldBuffer("rhoE",&rhoE[0]);

  simState2d.SetFieldBuffer("pressure",&p[0]);
  simState2d.SetFieldBuffer("temperature",&T[0]);
  simState2d.SetFieldBuffer("velocity",&V[0]);
  simState2d.SetFieldBuffer("rhom1",&rhom1[0]);

  state_t rhsState(simState2d);


  // -- Set up the RHS --
  rhs_t viscidRHS2d;
  
  if(viscidRHS2d.Initialize(simGrid2d,simState2d,paramState2d,sbpOperator2d)){
    parTestViscidRHS = false;
    std::cout << "Failed to initialize RHS in TestViscidRHS" << std::endl;
    return;
  }
 

  std::vector<int> numThreadsDim;
  simGrid2d.SetupThreads(numThreadsDim);
  viscidRHS2d.InitThreadIntervals();

  if(viscidRHS2d.SetupRHS(simGrid2d,simState2d,rhsState)){
    parTestViscidRHS = false;
    std::cout << "Failed to setup RHS in TestViscidRHS" << std::endl;
    return;
  }


  parallelUnitResults.UpdateResult("Viscid:RHS:RunsInParallel",parTestViscidRHS);

  // 
  // setup a 2D periodic velocity field so we can test the velocity gradient
  // setup a 2D periodic temperature field so we can test the temperature gradient
  // leave the density alone, update the energy
  //
  std::vector<double> velGradExpected(numDim*numDim*numPointsBuffer,-999.0);
  std::vector<double> tempGradExpected(numDim*numPointsBuffer,-999.0);
  std::vector<double> uAmp; 
  std::vector<double> vAmp;
  std::vector<double> tempAmp;
  uAmp.push_back(1.0);
  uAmp.push_back(2.0);
  uAmp.push_back(0.0);
  vAmp.push_back(3.0);
  vAmp.push_back(4.0);
  vAmp.push_back(0.0);
  tempAmp.push_back(5.0);
  tempAmp.push_back(6.0);
  tempAmp.push_back(0.0);
  std::vector<double> period(numDim,0);
  for(int i=0;i<numDim;i++){
    period[i]=dX[i]*numNodes;
  }
  double x0 = 0.;
  double y0 = 0.;
  size_t iStart = partitionBufferInterval2d[0].first;
  size_t iEnd   = partitionBufferInterval2d[0].second;
  size_t jStart = partitionBufferInterval2d[1].first;
  size_t jEnd   = partitionBufferInterval2d[1].second;
  for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
    size_t jBufferIndex = jIndex*bufferSizes[0];
    size_t jPartIndex = (jIndex - jStart) + partitionInterval2d[1].first;
    double y = jPartIndex*dX[1];
    for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
      size_t bufferIndex = jBufferIndex + iIndex;
      size_t iPartIndex = (iIndex - iStart) + partitionInterval2d[0].first;
      double x = iPartIndex*dX[0];

      rhoV[bufferIndex] = rho[bufferIndex]*
        testfixtures::viscid::Cosine(uAmp,period,x0,y0,0.,x,y,0.,-1);
      rhoV[bufferIndex+numPointsBuffer] = rho[bufferIndex]*
        testfixtures::viscid::Cosine(vAmp,period,x0,y0,0.,x,y,0.,-1);

      // just the internal part, add a linear offset
      rhoE[bufferIndex] += rho[bufferIndex]/(gamma-1)*
        testfixtures::viscid::Cosine(tempAmp,period,x0,y0,0.,x,y,0.,-1);

      //std::cout << "Expected velocity (" << rhoV[bufferIndex] << "," 
                //<< rhoV[bufferIndex+numPointsBuffer] << "," 
                //<< rhoV[bufferIndex+2*numPointsBuffer] 
                //<< ") at (" << x << "," << y << "," << z << ")" << std::endl;
      rhoE[bufferIndex] = rhoE[bufferIndex] + (rhoV[bufferIndex]*rhoV[bufferIndex] +
                             rhoV[bufferIndex+numPointsBuffer]*
                             rhoV[bufferIndex+numPointsBuffer])/(2.0*rho[bufferIndex]);

      // calculate the exact velocity gradient for comparison
      // one big flat vector with numPointsBuffer*numDin^2 size
      // each {} is numPoints in size
      // [ {du/dx} {dv/dx} {du/dy} {dv/dy} ]
      for(int dIndex = 0; dIndex<numDim; dIndex++){
        velGradExpected[bufferIndex+dIndex*numDim*numPointsBuffer] = 
          testfixtures::viscid::Cosine(uAmp,period,x0,y0,0.,x,y,0.,dIndex)*dX[dIndex];
        velGradExpected[bufferIndex+numPointsBuffer+dIndex*numDim*numPointsBuffer] = 
          testfixtures::viscid::Cosine(vAmp,period,x0,y0,0.,x,y,0.,dIndex)*dX[dIndex];
      }

      // calculate the exact temperature gradient for comparison
      // one big flat vector with numPointsBuffer*numDim size
      // each {} is numPoints in size
      // [ {dT/dx} {dT/dy} ]
      // we remove the grid space as the metrics are applied inside ComputeStressTensor
      for(int dIndex = 0; dIndex<numDim; dIndex++){
        tempGradExpected[bufferIndex+dIndex*numPointsBuffer] = 
          testfixtures::viscid::Cosine(tempAmp,period,x0,y0,0.,x,y,0.,dIndex)*dX[dIndex];
      }
    }
  }



  int returnCode = viscidRHS2d.ComputeDV(0);
  if(returnCode > 0){
    VelocityGradientWorks = false;
    std::cout << "ComputeVelGrad failed when computing DV " << std::endl;
    return;
  }
  //std::cout << "ComputeDV done" << std::endl;


  std::cout << "Calculating Velocity Gradient" << std::endl;
  returnCode = viscidRHS2d.ComputeVelGrad(0);
  if(returnCode > 0){
    std::cout << "ComputeVelGrad failed, aborting tests." << std::endl;
    VelocityGradientWorks = false;
    return;
  }

  simGrid2d.Communicator().Barrier();

  VelocityGradientWorks = CheckResult(VelocityGradientWorks,testComm);
  if(!VelocityGradientWorks){
    std::cout << "Velocity Gradient Calculation Failed" << std::endl;
  }

  // check the velocity gradient
  // first get the gradient from the rhs object
  // then check against the expected answer
  simGrid2d.Communicator().Barrier();

  if(VelocityGradientWorks){
    std::cout << "Checking Velocity Gradient calculation" << std::endl;
  
    const std::vector<double *> &velGradBuffers(viscidRHS2d.VelGradBuffers());
    
    bool printVelGrad = true;
    for(int iDim = 0;iDim < numDim;iDim++){
      const double* velGradData = velGradBuffers[iDim];
      
      for(int jDim = 0;jDim < numDim;jDim++){
        const double* velGradDataComponent = &velGradData[jDim*numPointsBuffer];
        
        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jBufferIndex = jIndex*bufferSizes[0];
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jBufferIndex + iIndex;
            size_t expectedIndex = numPointsBuffer*iDim*numDim+jDim*numPointsBuffer+bufferIndex;
            
            
            double error = std::abs(velGradDataComponent[bufferIndex] - 
                                    velGradExpected[expectedIndex])/velGradExpected[expectedIndex];
            
            if(error > 1e-2){
              VelocityGradientWorks = false;
              if(printVelGrad) {
                std::cout << std::scientific << std::setprecision(16)
                          << "In ComputeVelGrad: "<< std::endl
                          << "Expected velGrad=" << velGradExpected[expectedIndex] << std::endl
                          << " got velGrad=" << velGradDataComponent[bufferIndex] << std::endl
                          << " at (i,j,bufferIndex)=(" << iDim << "," << jDim << "," 
                          << bufferIndex << ")" << std::endl;
                std::cout << "Error =" << error << std::endl;;
                std::cout << "No further error messages will be generated for velGrad 2D Cosine" << std::endl;
                printVelGrad=false;
              }
            }
          }
        }
      }
    }
  }

  VelocityGradientWorks = CheckResult(VelocityGradientWorks,testComm);
  if(!VelocityGradientWorks){
    std::cout << "Velocity Gradient Compare Failed" << std::endl;
  }

  // Compute the temperature gradient
  std::cout << "Calculating Temperature Gradient" << std::endl;
  returnCode = viscidRHS2d.ComputeTemperatureGrad(0);
  if(returnCode > 0){
    TemperatureGradientWorks = false;
    return;
  }

  TemperatureGradientWorks = CheckResult(TemperatureGradientWorks,testComm);
  if(!TemperatureGradientWorks){
    std::cout << "Temperature Gradient Calculation Failed" << std::endl;
  }

  // check the temperature gradient
  // first get the gradient from the rhs object
  // then check against the expected answer
  std::cout << "Checking Temperature Gradient calculation" << std::endl;
  const std::vector<double *> &tempGradBuffers(viscidRHS2d.TemperatureGradBuffers());

  if(TemperatureGradientWorks){
    bool printTempGrad = true;
    for(int iDim = 0;iDim < numDim;iDim++){
      const double* tempGradDataComponent = tempGradBuffers[iDim];
      
      for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
        size_t jBufferIndex = jIndex*bufferSizes[0];
        for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
          size_t bufferIndex = jBufferIndex + iIndex;
          size_t expectedIndex = bufferIndex+numPointsBuffer*iDim;
          
          double error = std::abs(tempGradDataComponent[bufferIndex] - 
            tempGradExpected[expectedIndex])/tempGradExpected[expectedIndex];
          
          if(error > 1e-2){
            TemperatureGradientWorks = false;
            if(printTempGrad) {
              std::cout << std::scientific << std::setprecision(16)
                        << "In ComputeTemperatureGrad: "<< std::endl
                        << "Expected tempGrad=" << tempGradExpected[expectedIndex] << std::endl
                        << " got tempGrad=" << tempGradDataComponent[bufferIndex] << std::endl
                        << " at (i,bufferIndex)=(" << iDim << "," 
                        << bufferIndex << ")" << std::endl;
              std::cout << "Temperature " << T[bufferIndex] << std::endl;
              std::cout << "Error =" << error << std::endl;;
              std::cout << "No further error messages will be generated for tempGrad 2D periodic" 
                        << std::endl;
              printTempGrad=false;
            }
          }
        }
      }
    }
  }
  TemperatureGradientWorks = CheckResult(TemperatureGradientWorks,testComm);
  if(!TemperatureGradientWorks){
    std::cout << "Temperature Gradient Compare Failed" << std::endl;
  }

  std::cout << "Done with 2D periodic Velocity Gradient Test" << std::endl;
  std::cout << "Starting 3D periodic Velocity Gradient Test" << std::endl;

  //numNodes = 51;
  gridSizes.resize(3,numNodes);
  simGrid3d.SetGridSizes(gridSizes);

  physicalExtent.resize(6,1.0);
  physicalExtent[0] = 0.0;
  physicalExtent[2] = 0.0;
  physicalExtent[4] = 0.0;
  simGrid3d.SetPhysicalExtent(physicalExtent);

  periodicDirs.resize(3,true);
  simGrid3d.SetPeriodicDirs(periodicDirs);
   
  if(euler::util::InitializeSimulationFixtures(simGrid3d,sbpOperator3d,2,testComm,&messageStream)){
    parTestViscidRHS = false;
    std::cout << "Failed to setup 3D grid in TestViscidRHS" << std::endl;
    return;
  }

  std::cout << messageStream.str() << std::endl;
  messageStream.str("");
  messageStream.clear();

  const pcpp::IndexIntervalType &partitionBufferInterval3d(simGrid3d.PartitionBufferInterval());
  const pcpp::IndexIntervalType &partitionInterval3d(simGrid3d.PartitionInterval());
  numDim = gridSizes.size();

  if(testfixtures::viscid::SetupViscidState(simGrid3d,simState3d,paramState3d,0)){
    parTestViscidRHS = false;
    std::cout << "Failed to setup state in TestViscidRHS" << std::endl;
    return;
  }

  numPointsBuffer = simGrid3d.BufferSize();
  //bufferSizes.resize(simGrid3d.BufferSizes());
  const std::vector<size_t> &bufferSizes3d(simGrid3d.BufferSizes());

  //std::vector<double> &dX.resize(simGrid3d.DX());
  std::vector<double> &dX3d(simGrid3d.DX());

  // Add some fields that are not usual
  paramState3d.AddField("dX",'s',numDim,8,"space");
  paramState3d.Create(numPointsBuffer,0);

  paramState3d.SetFieldBuffer("gamma",&gamma);
  paramState3d.SetFieldBuffer("inputDT",&dt);
  paramState3d.SetFieldBuffer("inputCFL",&cfl);
  paramState3d.SetFieldBuffer("dX",&dX3d[0]);
  paramState3d.SetFieldBuffer("refRe",&Re);
  paramState3d.SetFieldBuffer("refPr",&Pr);
  paramState3d.SetFieldBuffer("power",&power);
  paramState3d.SetFieldBuffer("beta",&beta);
  paramState3d.SetFieldBuffer("bulkViscFac",&bulkViscFac);

  rho.resize(numPointsBuffer,1.0);
  rhoV.resize(numDim*numPointsBuffer,0.0);
  rhoE.resize(numPointsBuffer);
  rhoE.assign(rhoE.size(),1.0/(gamma-1));

  T.resize(numPointsBuffer,0.0);
  p.resize(numPointsBuffer,0.0);
  V.resize(numDim*numPointsBuffer,0.0);

  rhom1.resize(numPointsBuffer,0.0);

  simState3d.SetFieldBuffer("simTime",&t);
  simState3d.SetFieldBuffer("rho",&rho[0]);
  simState3d.SetFieldBuffer("rhoV",&rhoV[0]);
  simState3d.SetFieldBuffer("rhoE",&rhoE[0]);

  simState3d.SetFieldBuffer("pressure",&p[0]);
  simState3d.SetFieldBuffer("temperature",&T[0]);
  simState3d.SetFieldBuffer("velocity",&V[0]);
  simState3d.SetFieldBuffer("rhom1",&rhom1[0]);

  state_t rhsState3d(simState3d);

  // -- Set up the RHS --
  rhs_t viscidRHS3d;
  
  if(viscidRHS3d.Initialize(simGrid3d,simState3d,paramState3d,sbpOperator3d)){
    parTestViscidRHS = false;
    std::cout << "Failed to initialize RHS in TestViscidRHS" << std::endl;
    return;
  }

  std::vector<int> numThreadsDim3;
  simGrid3d.SetupThreads(numThreadsDim3);
  viscidRHS3d.InitThreadIntervals();

  if(viscidRHS3d.SetupRHS(simGrid3d,simState3d,rhsState3d)){
    parTestViscidRHS = false;
    std::cout << "Failed to setup RHS in TestViscidRHS" << std::endl;
    return;
  }

  parallelUnitResults.UpdateResult("Viscid:RHS:RunsInParallel",parTestViscidRHS);

  // 
  // setup a 3D periodic velocity field so we can test the velocity gradient
  // setup a 3D periodic temperature field so we can test the temperature gradient
  // leave the density alone, update the energy even though we don't really need it
  //
  velGradExpected.resize(numDim*numDim*numPointsBuffer,-999.0);
  tempGradExpected.resize(numDim*numPointsBuffer,-999.0);
  uAmp.resize(0); 
  vAmp.resize(0);
  tempAmp.resize(0);
  std::vector<double> wAmp;
  uAmp.push_back(1.0);
  uAmp.push_back(2.0);
  uAmp.push_back(3.0);
  vAmp.push_back(4.0);
  vAmp.push_back(5.0);
  vAmp.push_back(6.0);
  wAmp.push_back(7.0);
  wAmp.push_back(8.0);
  wAmp.push_back(9.0);
  tempAmp.push_back(10.0);
  tempAmp.push_back(11.0);
  tempAmp.push_back(12.0);
  period.resize(numDim,0);
  for(int i=0;i<numDim;i++){
    period[i]=dX3d[i]*numNodes;
  }
  x0 = 0.;
  y0 = 0.;
  double z0 = 0.;
  iStart = partitionBufferInterval3d[0].first;
  iEnd   = partitionBufferInterval3d[0].second;
  jStart = partitionBufferInterval3d[1].first;
  jEnd   = partitionBufferInterval3d[1].second;
  size_t kStart = partitionBufferInterval3d[2].first;
  size_t kEnd   = partitionBufferInterval3d[2].second;
  size_t nPlane = bufferSizes3d[0]*bufferSizes3d[1];
  for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
    size_t kBufferIndex = kIndex*nPlane;
    size_t kPartIndex = (kIndex - kStart) + partitionInterval3d[2].first;
    double z = kPartIndex*dX3d[2]; 
    for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
      size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes3d[0];
      size_t jPartIndex = (jIndex - jStart) + partitionInterval3d[1].first;
      double y = jPartIndex*dX3d[1];
      for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
        size_t bufferIndex = jkBufferIndex + iIndex;
        size_t iPartIndex = (iIndex - iStart) + partitionInterval3d[0].first;
        double x = iPartIndex*dX3d[0];

        rhoV[bufferIndex] = rho[bufferIndex]*
          testfixtures::viscid::Cosine(uAmp,period,x0,y0,z0,x,y,z,-1);
        rhoV[bufferIndex+numPointsBuffer] = rho[bufferIndex]*
          testfixtures::viscid::Cosine(vAmp,period,x0,y0,z0,x,y,z,-1);
        rhoV[bufferIndex+2*numPointsBuffer] = rho[bufferIndex]*
          testfixtures::viscid::Cosine(wAmp,period,x0,y0,z0,x,y,z,-1);

        // just the internal part, add a linear offset
        rhoE[bufferIndex] += rho[bufferIndex]/(gamma-1)*
          testfixtures::viscid::Cosine(tempAmp,period,x0,y0,z0,x,y,z,-1);

        rhoE[bufferIndex] = rhoE[bufferIndex] + (rhoV[bufferIndex]*rhoV[bufferIndex] +
                               rhoV[bufferIndex+numPointsBuffer]*
                               rhoV[bufferIndex+numPointsBuffer] +
                               rhoV[bufferIndex+2*numPointsBuffer]*
                               rhoV[bufferIndex+2*numPointsBuffer])/(2.0*rho[bufferIndex]);

        // calculate the exact velocity gradient for comparison
        // one big flat vector with numPointsBuffer*numDin^2 size
        // each {} is numPoints in size
        // [ {du/dx} {dv/dx} {dw/dx} {du/dy} {dv/dy} {dw/dy} {du/dz} {dv/dz} {dw/dz} ]
        for(int dIndex = 0; dIndex<numDim; dIndex++){
          velGradExpected[bufferIndex+dIndex*numDim*numPointsBuffer] = 
            testfixtures::viscid::Cosine(uAmp,period,x0,y0,z0,x,y,z,dIndex)*dX3d[dIndex];
          velGradExpected[bufferIndex+numPointsBuffer+dIndex*numDim*numPointsBuffer] = 
            testfixtures::viscid::Cosine(vAmp,period,x0,y0,z0,x,y,z,dIndex)*dX3d[dIndex];
          velGradExpected[bufferIndex+2*numPointsBuffer+dIndex*numDim*numPointsBuffer] = 
            testfixtures::viscid::Cosine(wAmp,period,x0,y0,z0,x,y,z,dIndex)*dX3d[dIndex];
        }

        // calculate the exact temperature gradient for comparison
        // one big flat vector with numPointsBuffer*numDim size
        // each {} is numPoints in size
        // [ {dT/dx} {dT/dy} {dT/dz}]
        // we remove the grid space as the metrics are applied inside ComputeStressTensor
        for(int dIndex = 0; dIndex<numDim; dIndex++){
          tempGradExpected[bufferIndex+dIndex*numPointsBuffer] = 
            testfixtures::viscid::Cosine(tempAmp,period,x0,y0,z0,x,y,z,dIndex)*dX3d[dIndex];
        }
      }
    }
  }

  returnCode = viscidRHS3d.ComputeDV(0);
  if(returnCode > 0){
    VelocityGradientWorks = false;
    std::cout << "ComputeVelGrad failed when computing DV " << std::endl;
    return;
  }

  std::cout << "Calling ComputeVelGrad" << std::endl;
  returnCode = viscidRHS3d.ComputeVelGrad(0);
  if(returnCode > 0){
    VelocityGradientWorks = false;
    return;
  }

  VelocityGradientWorks = CheckResult(VelocityGradientWorks,testComm);
  if(!VelocityGradientWorks){
    std::cout << "Velocity Gradient Calculation Failed" << std::endl;
  }

  if(VelocityGradientWorks){
    std::cout << "Checking Velocity Gradient calculation" << std::endl;
    // check the velocity gradient
    // first get the gradient from the rhs object
    // then check against the expected answer
    const std::vector<double *> &velGradBuffers3d(viscidRHS3d.VelGradBuffers());
    
    bool printVelGrad = true;
    for(int iDim = 0;iDim < numDim;iDim++){
      const double* velGradData3d = velGradBuffers3d[iDim];
      
      for(int jDim = 0;jDim < numDim;jDim++){
        const double* velGradDataComponent3d = &velGradData3d[jDim*numPointsBuffer];
        
        for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
          size_t kBufferIndex = kIndex*nPlane;
          for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
            size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes3d[0];
            for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
              size_t bufferIndex = jkBufferIndex + iIndex;
              size_t expectedIndex = numPointsBuffer*iDim*numDim+jDim*numPointsBuffer+bufferIndex;
              
              double error = std::abs(velGradDataComponent3d[bufferIndex] - 
                                      velGradExpected[expectedIndex])/velGradExpected[expectedIndex];
              
              if(error > 1e-2){
                VelocityGradientWorks = false;
                if(printVelGrad) {
                  std::cout << std::scientific << std::setprecision(16)
                            << "In ComputeVelGrad: "<< std::endl
                            << "Expected velGrad=" << velGradExpected[expectedIndex] << std::endl
                            << " got velGrad=" << velGradDataComponent3d[bufferIndex] << std::endl
                            << " at (i,j,bufferIndex)=(" << iDim << "," << jDim << "," 
                            << bufferIndex << ")" << std::endl;
                  std::cout << "Error =" << error << std::endl;;
                  std::cout << "No further error messages will be generated for velGrad 3d" << std::endl;
                  printVelGrad=false;
                }
              }
            }
          }
        }
      }
    }
  }

  VelocityGradientWorks = CheckResult(VelocityGradientWorks,testComm);
  if(!VelocityGradientWorks){
    std::cout << "Velocity Gradient Compare Failed" << std::endl;
  }

  // Compute the temperature gradient
  std::cout << "Calculating Temperature Gradient" << std::endl;
  returnCode = viscidRHS3d.ComputeTemperatureGrad(0);
  if(returnCode > 0){
    TemperatureGradientWorks = false;
    return;
  }

  TemperatureGradientWorks = CheckResult(TemperatureGradientWorks,testComm);
  if(!TemperatureGradientWorks){
    std::cout << "Temperature Gradient Calculation Failed" << std::endl;
  }

  // check the temperature gradient
  // first get the gradient from the rhs object
  // then check against the expected answer
  std::cout << "Checking Temperature Gradient calculation" << std::endl;
  const std::vector<double *> &tempGradBuffers3d(viscidRHS3d.TemperatureGradBuffers());

  if(TemperatureGradientWorks){
    bool printTempGrad = true;
    for(int iDim = 0;iDim < numDim;iDim++){
      const double* tempGradDataComponent = tempGradBuffers3d[iDim];

      for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
        size_t kBufferIndex = kIndex*nPlane;
        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes3d[0];
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jkBufferIndex + iIndex;
            size_t expectedIndex = bufferIndex + numPointsBuffer*iDim;
    
            double error = std::abs(tempGradDataComponent[bufferIndex] - 
              tempGradExpected[expectedIndex])/tempGradExpected[expectedIndex];
          
            if(error > 1e-2){
              TemperatureGradientWorks = false;
              if(printTempGrad) {
                std::cout << std::scientific << std::setprecision(16)
                          << "In ComputeTemperatureGrad: "<< std::endl
                          << "Expected tempGrad=" << tempGradExpected[expectedIndex] << std::endl
                          << " got tempGrad=" << tempGradDataComponent[bufferIndex] << std::endl
                          << " at (i,bufferIndex)=(" << iDim << "," 
                          << bufferIndex << ")" << std::endl;
                std::cout << "Temperature " << T[bufferIndex] << std::endl;
                std::cout << "Error =" << error << std::endl;;
                std::cout << "No further error messages will be generated for tempGrad 3D Linear" 
                          << std::endl;
                printTempGrad=false;
              }
            }
          }
        }
      }
    }
  }

  std::cout << "Done with 3D periodic Velocity Gradient Test" << std::endl;
  
  parallelUnitResults.UpdateResult("Viscid:RHS:VelocityGradientPeriodic",VelocityGradientWorks);
  parallelUnitResults.UpdateResult("Viscid:RHS:TemperatureGradientPeriodic",TemperatureGradientWorks);
  
  std::cout << "Done TestVelocityGradientPeriodic..." << std::endl;
  
};

void TestViscidRHS(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{
  std::cout << "TestViscidRHS..." << std::endl;

  typedef plascom2::operators::sbp::base  operator_t;
  typedef simulation::grid::halo          halo_t;
  typedef simulation::state::base         state_t;
  typedef pcpp::ParallelGlobalType        global_t;
  typedef pcpp::CommunicatorType          comm_t;
  typedef pcpp::IndexIntervalType         interval_t;

  typedef simulation::grid::parallel_blockstructured         grid_t;
  typedef simulation::domain::base<grid_t,state_t>  domain_t;
  typedef euler::rhs<grid_t,state_t,operator_t>             rhs_t;

  std::ostringstream messageStream;

  grid_t     simGrid;
  operator_t sbpOperator;
  halo_t     &simHalo(simGrid.Halo());
  state_t    simState;
  state_t    paramState;
  
  bool parTestViscidRHS = true;
  std::vector<size_t> gridSizes(3,4);

  if(testfixtures::CreateSimulationFixtures(simGrid,simHalo,sbpOperator,
                                            gridSizes,2,testComm,&messageStream)){
    parTestViscidRHS = false;
    return;
  }

  pcpp::IndexIntervalType &partitionBufferInterval(simGrid.PartitionBufferInterval());

  if(testfixtures::viscid::SetupViscidState(simGrid,simState,paramState,0)){
    parTestViscidRHS = false;
    return;
  }
  
  std::cout << messageStream.str() << std::endl;
  
  size_t numPointsBuffer = simGrid.BufferSize();

  double t = 0.0;
  double gamma = 1.4;
  double power = 2.0/3.0;
  double beta = 1.0;
  double bulkViscFac = 0.6;
  double dt = .001;
  double cfl = 1.0;
  double Re  = 100.0;
  double Pr  = 0.75;
  //  double numbers[] = {};
  int    flags = 1;

  int numDim = 3;

  paramState.SetFieldBuffer("gamma",&gamma);
  paramState.SetFieldBuffer("inputDT",&dt);
  paramState.SetFieldBuffer("inputCFL",&cfl);
  paramState.SetFieldBuffer("refRe",&Re);
  paramState.SetFieldBuffer("refPr",&Pr);
  paramState.SetFieldBuffer("power",&power);
  paramState.SetFieldBuffer("beta",&beta);
  paramState.SetFieldBuffer("bulkViscFac",&bulkViscFac);
  //  paramState.SetFieldBuffer("Numbers",numbers);
  paramState.SetFieldBuffer("Flags",&flags);

  std::vector<double> rho(numPointsBuffer,1.0);
  std::vector<double> rhoV(3*numPointsBuffer,0.0);
  std::vector<double> rhoE(3*numPointsBuffer,1.0/(gamma-1.0));

  std::vector<double> dRho(numPointsBuffer,0.0);
  std::vector<double> dRhoV(3*numPointsBuffer,0.0);
  std::vector<double> dRhoE(numPointsBuffer,0.0);

  std::vector<double> T(numPointsBuffer,0.0);
  std::vector<double> p(numPointsBuffer,0.0);
  std::vector<double> V(3*numPointsBuffer,0.0);

  std::vector<double> rhom1(numPointsBuffer,0.0);


  simState.SetFieldBuffer("simTime",&t);
  simState.SetFieldBuffer("rho",&rho[0]);
  simState.SetFieldBuffer("rhoV",&rhoV[0]);
  simState.SetFieldBuffer("rhoE",&rhoE[0]);

  simState.SetFieldBuffer("pressure",&p[0]);
  simState.SetFieldBuffer("temperature",&T[0]);
  simState.SetFieldBuffer("velocity",&V[0]);
  simState.SetFieldBuffer("rhom1",&rhom1[0]);

  state_t rhsState(simState);

  rhsState.SetFieldBuffer("rho",&dRho[0]);
  rhsState.SetFieldBuffer("rhoV",&dRhoV[0]);
  rhsState.SetFieldBuffer("rhoE",&dRhoE[0]);

  // set up the domain - don't care about this in testing
  //   domain_t eulerDomain;
  //   eulerDomain.numGrids = 1;
  //   eulerDomain.Grids().resize(1,&simGrid);
  //   eulerDomain.States().resize(1,&simState);

  // -- Set up the RHS --
  rhs_t eulerRHS;
  
  std::cout << "Initialize" << std::endl;
  if(eulerRHS.Initialize(simGrid,simState,paramState,sbpOperator)){
    parTestViscidRHS = false;
    return;
  }

  std::cout << "SetupRHS" << std::endl;
  if(eulerRHS.SetupRHS(simGrid,simState,rhsState)){
    parTestViscidRHS = false;
    return;
  }

  simGrid.SetNumThreads(1);
  eulerRHS.InitThreadIntervals();
 
  
  std::cout << "RHS" << std::endl;
  int returnCode = eulerRHS.RHS(0);
  if(returnCode > 0){
    parTestViscidRHS = false;
  }
  std::cout << "Done RHS" << std::endl;

  parTestViscidRHS = CheckResult(parTestViscidRHS,testComm);
  if(!parTestViscidRHS){
    std::cout << "Problem 0" << std::endl;
  }

  if(parTestViscidRHS){
    for(size_t iPoint = 0;(iPoint < numPointsBuffer)&&parTestViscidRHS;iPoint++){
      if(std::abs(dRho[iPoint]) > 1e-15) parTestViscidRHS = false;
      for(int iDim = 0;iDim < numDim;iDim++)
        if(std::abs(dRhoV[iPoint+iDim*numPointsBuffer]) > 1e-15) parTestViscidRHS = false;
      if(std::abs(dRhoE[iPoint]) > 1e-15) parTestViscidRHS = false;
    }
  }

  parTestViscidRHS = CheckResult(parTestViscidRHS,testComm);
  if(!parTestViscidRHS){
    std::cout << "Problem 1" << std::endl;
  }

  parallelUnitResults.UpdateResult("Viscid:ViscidRHS:3dWorksInParallel",parTestViscidRHS);


};

void TestViscidKernelsCurvilinear(ix::test::results &serialUnitResults)
{
  
  std::cout << "Running viscid kernels tests with full curvilinear metrics." 
            << std::endl;

  bool tauCurvilinear = true;
  bool qCurvilinear   = true;
  bool flux1dTest     = true;

  for(int numDim = 1;numDim <= 3;numDim++){


    std::vector<size_t> numPoints(numDim,1);
    size_t numPointsBuffer = 1;
    for(int iDim = 0;iDim < numDim;iDim++){
      numPoints[iDim] = 2*std::pow(2,iDim) + 1;
      numPointsBuffer *= numPoints[iDim];
      //      std::cout << "NumPoints[" << iDim << "] = " << numPoints[iDim] << std::endl;
    }

    //    std::cout << "NumPoints = " << numPointsBuffer << std::endl;

    std::vector<double> mu(numPointsBuffer,0.0);
    std::vector<double> lambda(numPointsBuffer,0.0);
    std::vector<double> gradVIJK(numDim*numDim*numPointsBuffer,0.0);
    std::vector<double> gradVXYZ(numDim*numDim*numPointsBuffer,0.0);
    std::vector<double> gridJacobian(numPointsBuffer,0.0);
    std::vector<double> gridMetric(numDim*numDim*numPointsBuffer,0.0);
    std::vector<double> velocity(numDim*numPointsBuffer,0.0);
    std::vector<double> gradT(numDim*numPointsBuffer,0.0);
    std::vector<double> tau(numDim*numDim*numPointsBuffer,0.0);
    std::vector<double> expectedTau(numDim*numDim*numPointsBuffer,0.0);

    double piValue = 4.0*std::atan(1.0);
    double Re = 25.0;
    std::vector<double> metricValues(numDim*numDim,0.0);
    std::vector<double> gradValuesIJK(numDim*numDim,0.0);
    std::vector<double> gradValuesXYZ(numDim*numDim,0.0);
    // All the metric and gradient components are set to 
    // the following for all points:
    //
    // metric =   |.01 .02 |   | .01  .02 .03 |
    //            |.03 .04 |,  | .04  .05 .06 |
    //                         | .07  .08 .09 |
    //          
    // **NOTE: gradVIJK goes like: dV(j)/dX(i) (i.e. the *transpose* of gradV(ijk))
    //
    // gradVIJK  =   | 1  2 |     | 1 2 3 |
    //               | 3  4 |,    | 4 5 6 |
    //                            | 7 8 9 |
    for(int iDim = 0;iDim < numDim;iDim++){
      for(int jDim = 0;jDim < numDim;jDim++){
        int component = iDim*numDim + jDim;
        gradValuesIJK[component] = static_cast<double>(component + 1);
        metricValues[component] = static_cast<double>(component+1)/100.0;
        //         std::cout << "GradValuesIJK[" << iDim << "][" << jDim << "] = " << gradValuesIJK[component] << std::endl
        //                   << "metricValues[" << iDim << "][" << jDim << "] = " << metricValues[component] << std::endl;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
          gradVIJK[component*numPointsBuffer+iPoint] = gradValuesIJK[component];
          gridMetric[component*numPointsBuffer+iPoint] = metricValues[component];
        }
      }
    }
    
    // Velocity = <1, 2, 3>
    // gradT    = <4 , 5, 6>
    for(int iDim = 0;iDim < numDim;iDim++){
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        velocity[iDim*numPointsBuffer+iPoint] = static_cast<double>(iDim+1);
        gradT[iDim*numPointsBuffer+iPoint] = static_cast<double>(iDim+4);
      }
    }
    
    // Jacobian = iPoint + 1 
    // mu = ((iPoint*iPoint) + 1)/1000.0
    // lambda = COS(iPoint/2PI)
    for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
      gridJacobian[iPoint] = static_cast<double>(iPoint+1);
      mu[iPoint] = static_cast<double>(iPoint*iPoint + 1)/1000.0;
      lambda[iPoint] = std::cos(static_cast<double>(iPoint+1)/(2.0*piValue));
    }

    
    
    double divValue = 0.0;
    for(int iDim = 0;iDim < numDim;iDim++){
      int gradOffset = iDim;
      for(int jDim = 0;jDim < numDim;jDim++){
        int component = iDim*numDim + jDim;
        int metricOffset = jDim;
        for(int kDim = 0;kDim < numDim;kDim++){
          gradValuesXYZ[component] += gradValuesIJK[gradOffset+kDim*numDim]*metricValues[jDim+kDim*numDim];
        }
        if(iDim == jDim)
          divValue += gradValuesXYZ[component];
      }
    }

    pcpp::IndexIntervalType regionInterval;
    regionInterval.InitSimple(numPoints);
    std::vector<size_t> fortranInterval;
    regionInterval.Flatten(fortranInterval);
    for(int iDim = 0;iDim < numDim;iDim++){
      fortranInterval[iDim*2]++;
      fortranInterval[iDim*2+1]++;
    }
    
    int gridType = simulation::grid::CURVILINEAR;
    std::vector<double*> bufferPtrArray(10,NULL);
    bufferPtrArray[0] = &mu[0];
    bufferPtrArray[1] = &lambda[0];
    bufferPtrArray[2] = &lambda[0];
    std::vector<double *> velGradBuffers(numDim,NULL);
    std::vector<double *> tauBuffers(numDim*numDim,NULL);
    
    for(int iDim = 0;iDim < numDim;iDim++){
      velGradBuffers[iDim] = &gradVIJK[iDim*numDim*numPointsBuffer];
      for(int jDim = 0;jDim < numDim;jDim++){
        size_t tauComponent = iDim*numDim + jDim;
        size_t tauTComponent = jDim*numDim + iDim;
        if(iDim <= jDim){
          tauBuffers[tauComponent] = &tau[tauComponent*numPointsBuffer];
          if(iDim == jDim){
            for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
              expectedTau[tauComponent*numPointsBuffer+iPoint] = (gridJacobian[iPoint]/Re)*
                (mu[iPoint]*(gradValuesXYZ[tauComponent]+gradValuesXYZ[tauTComponent]) +
                 lambda[iPoint]*divValue);
              //               std::cout << "expectedTau[" << iDim << "][" << jDim << "]["
              //                         << iPoint << "] = " << expectedTau[tauComponent*numPointsBuffer+iPoint]
              //                        << std::endl;
            }
          } else {
            for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
              expectedTau[tauComponent*numPointsBuffer+iPoint] = (gridJacobian[iPoint]/Re)*
                (mu[iPoint]*(gradValuesXYZ[tauComponent]+gradValuesXYZ[tauTComponent]));
              //               std::cout << "expectedTau[" << iDim << "][" << jDim << "]["
              //                         << iPoint << "] = " << expectedTau[tauComponent*numPointsBuffer+iPoint]
              //                         << std::endl;
            }
          }
        } else {
          tauBuffers[tauComponent] = &tau[tauTComponent*numPointsBuffer];
          for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
            expectedTau[tauComponent*numPointsBuffer+iPoint] = (gridJacobian[iPoint]/Re)*
              (mu[iPoint]*(gradValuesXYZ[tauComponent]+gradValuesXYZ[tauTComponent]));
            //             std::cout << "expectedTau[" << iDim << "][" << jDim << "]["
            //                       << iPoint << "] = " << expectedTau[tauComponent*numPointsBuffer+iPoint]
            //                       << std::endl;
            
          }
        }
      }
    }
    
    std::cout << "Checking tau in " << numDim << " dimensions...";
    if(viscid::util::ComputeTauBuffer(regionInterval,numPoints,gridType,gridMetric,gridJacobian,
                                      Re,bufferPtrArray,velGradBuffers,tauBuffers)){
      std::cout << "Compute Tau failed." << std::endl;
      return;
    }
    
    double errMax = -10000.0;
    double errMin = 1000000.0;
    double error = 0;
    for(int iDim = 0;iDim < numDim;iDim++){
      for(int jDim = 0;jDim < numDim;jDim++){
        size_t tauComponent = iDim*numDim + jDim;
        size_t tauIndex = tauComponent*numPointsBuffer;
        if(iDim <= jDim){
          for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
            error = std::abs(tau[tauComponent+iPoint] - 
                             expectedTau[tauComponent+iPoint])/std::abs(expectedTau[tauComponent+iPoint]);
            //             std::cout << "Actual Tau[" << iDim << "][" << jDim << "][" << iPoint << "] = "
            //                       << tau[tauComponent+iPoint] << std::endl;
            if(error > errMax) errMax = error;
            if(error < errMin) errMin = error;
          }
        }
      }
    }
    double tol = 1e-15;
    if(errMax > tol){
      tauCurvilinear = false;
      std::cout << "failed." << std::endl
                << "TAU Error: " << errMax << std::endl;
    } else {
      std::cout << "passed." << std::endl;
    }
    
    std::vector<double> expectedQ(numDim*numPointsBuffer,0.0);
    std::vector<double> q(numDim*numPointsBuffer,0.0);
    std::vector<double *> gradTBuffers(numDim,NULL);
    std::vector<double *> heatFluxBuffers(numDim,NULL);
    for(int iDim = 0;iDim < numDim;iDim++){
      gradTBuffers[iDim] = &gradT[iDim*numPointsBuffer];
      heatFluxBuffers[iDim] = &q[iDim*numPointsBuffer];
      size_t qIndex = iDim*numPointsBuffer;
      double dTdX = 0.0;
      for(int jDim = 0;jDim < numDim;jDim++){
        // Compute dT/dX(iDim)
        double gradTValue = gradT[jDim*numPointsBuffer];
        dTdX += gradTValue*metricValues[iDim+jDim*numDim];
      }

      for(int iPoint = 0;iPoint < numPointsBuffer;iPoint++)
        expectedQ[qIndex+iPoint] = lambda[iPoint]*gridJacobian[iPoint]*dTdX/Re;
    }

    std::cout << "Checking heat flux in " << numDim << " dimensions....";

    if(viscid::util::ComputeHeatFluxBuffer(regionInterval,numPoints,gridType,gridMetric,
                                           gridJacobian,Re,bufferPtrArray,gradTBuffers,
                                           heatFluxBuffers)){
      std::cout << "ComputeHeatFlux failed." << std::endl;
      return;
    }
    errMax = -100000.0;
    errMin = 1000000.0;
    for(int iPoint = 0;iPoint < numPointsBuffer;iPoint++){
      double err = std::abs(q[iPoint] - expectedQ[iPoint])/std::abs(expectedQ[iPoint]);
      if(err > tol){
        //        std::cout << "Q[" << iPoint << "] = " << q[iPoint] << " Expected: " << expectedQ[iPoint] << std::endl;
        qCurvilinear = false;
      }
      if(errMax < err) errMax = err;
      if(errMin > err) errMin = err;
    }
    if(errMax > tol){
      std::cout << "failed." << std::endl
                << "Max Heat Flux Error : " << errMax << std::endl;
    } else {
      std::cout << "passed." << std::endl;
    }

    // prescribe a whole other tau, one with whole numbers for easy debugging
    int tComponent = 0;
    int numTauComponents = (numDim*(numDim+1))/2;
    std::vector<double> tau2(numTauComponents*numPointsBuffer,0);
    for(int iDim = 0;iDim < numDim;iDim++){
      for(int jDim = 0;jDim < numDim;jDim++){
        if(iDim <= jDim){
          for(int iPoint = 0;iPoint < numPointsBuffer;iPoint++)
            tau2[tComponent*numPointsBuffer+iPoint] = mu[iPoint]*(tComponent+1);
          tComponent++;
        }
      }
    }

    std::vector<double> energyVector(numDim*numPointsBuffer,0);
    for(int iDim = 0;iDim < numDim;iDim++){
      for(int iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        energyVector[iDim*numPointsBuffer+iPoint] = lambda[iPoint]*(iDim+1);
      }
    }

    
    //    index2D = RESHAPE((/ 0, 1, 1, 2 /), SHAPE(index2D))
    //    index3D = RESHAPE((/ 0, 1, 2, 1, 3, 4, 2, 4, 5 /), SHAPE(index3D))
    std::vector<int> symmetricIndex(numDim*numDim,0);
    if(numDim == 1)
      symmetricIndex[0] = 0;
    if(numDim == 2){
      symmetricIndex[0] = 0;
      symmetricIndex[1] = 1;
      symmetricIndex[2] = 1;
      symmetricIndex[3] = 2;
    }
    if(numDim == 3){
      symmetricIndex[0] = 0;
      symmetricIndex[1] = 1;
      symmetricIndex[2] = 2;
      symmetricIndex[3] = 1;
      symmetricIndex[4] = 3;
      symmetricIndex[5] = 4;
      symmetricIndex[6] = 2;
      symmetricIndex[7] = 4;
      symmetricIndex[8] = 5;
    }

    std::cout << "Checking flux in " << numDim << " dimensions....";
    for(int iDim = 0;iDim < numDim;iDim++){
      
      int fluxDir = iDim+1;

      std::cout << "(" << fluxDir << ")...";
 
      int metricOffset = iDim*numDim*numPointsBuffer; // offset to appropriate row of metric

      std::vector<double> flux((numDim+2)*numPointsBuffer,0);
      std::vector<double> expectedFlux((numDim+2)*numPointsBuffer,0);

      FC_MODULE(viscid,strongflux1d,VISCID,STRONGFLUX1D)
        (&numDim,&fluxDir,&numPoints[0],&numPointsBuffer,&fortranInterval[0],
         &gridType,&gridMetric[0],&tau2[0],&energyVector[0],&flux[0]);
      
      // expected flux is row of metric dot column of tau (for momentum)
      for(int jDim = 0;jDim < numDim;jDim++){
        // need jDim(th) column of tau
        for(int kDim = 0;kDim < numDim;kDim++){
          int tauComponent = symmetricIndex[jDim+kDim*numDim];
          double *metricComponent = &gridMetric[metricOffset+kDim*numPointsBuffer];
          for(int iPoint = 0;iPoint < numPointsBuffer;iPoint++){
             expectedFlux[(jDim+1)*numPointsBuffer + iPoint] += 
               metricComponent[iPoint]*tau2[tauComponent*numPointsBuffer+iPoint];
          }
        }
      }

      // expected flux is row of metric dot Q (for energy)
      for(int jDim = 0;jDim < numDim;jDim++){
        double *metricComponent = &gridMetric[metricOffset+jDim*numPointsBuffer];
        for(int iPoint = 0;iPoint < numPointsBuffer;iPoint++){
          expectedFlux[(numDim+1)*numPointsBuffer + iPoint] +=
            metricComponent[iPoint]*energyVector[jDim*numPointsBuffer+iPoint];
        }
      }
      
      double maxError = -10000.0;
      double minError =  10000.0;
      for(int iDim = 0;iDim < numDim+2;iDim++){
        bool triggered = false;
        for(int iPoint = 0;iPoint < numPointsBuffer;iPoint++){
          double error = std::abs(flux[iDim*numPointsBuffer+iPoint] -
                                  expectedFlux[iDim*numPointsBuffer+iPoint])/
            std::abs(expectedFlux[iDim*numPointsBuffer+iPoint]);
          if(error > tol){
            if(!triggered){
              std::cout << "Flux component " << iDim << " failed." << std::endl;
              triggered = true;
            }
            flux1dTest = false;
          }
          if(error > maxError) maxError = error;
          if(error < minError) minError = error;
        }
      }
      if(maxError > tol){
        std::cout << "failed." << std::endl
                  << "Max flux error: " << maxError << std::endl;
        flux1dTest = false;
      }
    }
    if(flux1dTest){
      std::cout << "passed." << std::endl;
    }
  }

  serialUnitResults.UpdateResult("Viscid:ComputeHeatFlux:Curvilinear",qCurvilinear);
  serialUnitResults.UpdateResult("Viscid:ComputeTau:Curvilinear",tauCurvilinear);
  serialUnitResults.UpdateResult("Viscid:StrongFlux1D:Curvilinear",flux1dTest);

  return;

};