TestInviscid.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 "EulerTestFixtures.H"
#include "EOS.H"

using pcpp::comm::CheckResult;

typedef simulation::state::base state_t;

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

  typedef simulation::grid::parallel_blockstructured grid_t;

  std::cout << "TestEulerKernels.." << std::endl;
  
  // -- Set up the grid --

  //  int numDim = 2;

  std::vector<bool> dimComputeDV(2,true);
  std::vector<bool> dimUniformGridMetric(2,true);
  std::vector<bool> dimEulerRHS(2,true);
  std::vector<bool> dimUniformFlux(2,true);

  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);

    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];
    }
    testGrid.SetGridSpacings(dX);
    if(testGrid.Finalize()){
      std::cout << "Grid finalization failed." << std::endl;
      return;
    }
    // testGrid.DX() = dX[0];
    // testGrid.DY() = dX[1];
    
    // if(numDim == 3){
    //   testGrid.DZ() = dX[2];
    // }
    
    // testGrid.GlobalExtent().InitSimple(numNodes);
    // testGrid.LocalExtent().InitSimple(numNodes);
    // testGrid.NumberOfGlobalNodes() = numNodes;
    std::vector<size_t> bufferInterval;
    testGrid.PartitionBufferInterval().Flatten(bufferInterval);

    // std::cout << "Partition Buffer Interval: ";
    // testGrid.PartitionBufferInterval().PrettyPrint(std::cout);
    // std::cout << std::endl;

    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;
    double v = vVec[0];
    

    std::vector<double> expectedRhoM1(numNodesDomain,1.0);
    std::vector<double> expectedFlux((numDim+2)*numNodesDomain,0.0);
    std::vector<double> expectedPressure(numNodesDomain,1.0);
    std::vector<double> expectedTemperature(numNodesDomain,1.0);
    std::vector<double> expectedVHat(numDim*numNodesDomain,0.0);
    std::vector<double> expectedVelocity(numDim*numNodesDomain,0.0);
    std::vector<double> expectedInternalEnergy(numNodesDomain,1.0/(gamma-1.0));
    std::vector<double> expectedRHS((numDim+2)*numNodesDomain,0.0);
    
    std::vector<double> ke(numNodesDomain,0.0);
    std::vector<double> rho(numNodesDomain,1.0);
    std::vector<double> rhoV(numDim*numNodesDomain,1.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];
    }

    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];

    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 computeDVWorks = true;
    bool pfailed = false;
    bool tfailed = false;
    bool rfailed = false;
    bool vfailed = false;
    bool efailed = false;

    for(int iNode = 0;iNode < numNodesDomain;iNode++){
      if(std::abs(p[iNode] - expectedPressure[iNode])    > 1e-15){
        computeDVWorks = false;
        pfailed = true;
      }
      if(std::abs(T[iNode] - expectedTemperature[iNode]) > 1e-15){
        computeDVWorks = false;
        tfailed = true;
      }
      if(std::abs(rhom1[iNode] - expectedRhoM1[iNode])   > 1e-15){
        computeDVWorks = false;
        rfailed = true;
      }
      if(std::abs(internalEnergy[iNode] - expectedInternalEnergy[iNode])   > 1e-15){
        computeDVWorks = false;
        efailed = true;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        if(std::abs(velocity[iDim*numNodesDomain+iNode] - 
                    expectedVelocity[iDim*numNodesDomain+iNode]) > 1e-15){
          computeDVWorks = false;
          vfailed = true;
        }
      }
    }

    dimComputeDV[numDim-2] = computeDVWorks ;
    if(pfailed){
      std::cout << "Pressure failed" << std::endl;
    }
    if(tfailed){
      std::cout << "Temp failed" << std::endl;
    }
    if(rfailed){
      std::cout << "RhoM1 failed" << std::endl;
    }
    if(vfailed){
      std::cout << "Velo failed" << std::endl;
    }
    if(efailed){
      std::cout << "Internal energy failed" << std::endl;
    }


    // velocity.resize(0,0);
    // velocity.resize(numDim*numNodesDomain,1.0);
    std::vector<double> velHat(numDim*numNodesDomain,0.0);
    
    std::vector<size_t> fortranBufferInterval(2*numDim,1);
    for(int iDim = 0;iDim < numDim;iDim++){
      fortranBufferInterval[2*iDim+1] = bufferInterval[2*iDim+1] + 1;
    }
    
    FC_MODULE(grid,applyuniformgridmetric,GRID,APPLYUNIFORMGRIDMETRIC)
      (&numDim,&numNodes[0],&numNodesDomain,
       &fortranBufferInterval[0],
       &gridMetric[0],&velocity[0],&velHat[0]);

    std::vector<double *> vHatBuffers(numDim,(double *)NULL);
    for(int iDim = 0;iDim < numDim;iDim++){
      vBuffers[iDim] = &velocity[iDim*numNodesDomain];
      vHatBuffers[iDim] = &velHat[iDim*numNodesDomain];
    }
    //   double *vHat1 = &velHat[0];
    //   double *vHat2 = &velHat[numNodesDomain];
    //   double *vHat3 = &velHat[numNodesDomain*2];
    
    bool uniformGridMetricWorks = true;
    for(int iDim = 0;iDim < numDim;iDim++){
      size_t iOffset = iDim*numNodesDomain;
      for(int iNode = 0;iNode < numNodesDomain;iNode++){
        double error = std::abs(velHat[iOffset+iNode] - expectedVHat[iOffset+iNode]);
        if(error > 1e-15){
          uniformGridMetricWorks = false;
        }
      }
    }
    dimUniformGridMetric[numDim-2] = uniformGridMetricWorks;

    {
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t vectorOffset = iDim*numNodesDomain;
        std::vector<double> scaledPressure(numNodesDomain,0.0);
        for(size_t iPoint = 0;iPoint < numNodesDomain;iPoint++)
          scaledPressure[iPoint] = p[iPoint]*gridMetric[iDim];
        std::vector<double> eulerFlux(numNodesDomain*(numDim+2),0.0);
        std::vector<double> expectedFlux((numDim+2)*numNodesDomain,0.0);
        for(size_t iPoint = 0;iPoint < numNodesDomain;iPoint++){
          //expectedFlux[iPoint] = rho[iPoint]*velocity[vectorOffset+iPoint];
          expectedFlux[iPoint] = rho[iPoint]*velHat[vectorOffset+iPoint];
          size_t fOffset = 0;
          for(int ivDim = 0;ivDim < numDim;ivDim++){
            size_t iOffset = ivDim*numNodesDomain;
            fOffset = iOffset + numNodesDomain;
            expectedFlux[iPoint+fOffset] = 
              //rhoV[iOffset+iPoint]*velocity[vectorOffset+iPoint];
              rhoV[iOffset+iPoint]*velHat[vectorOffset+iPoint];
            if(ivDim == iDim)
              //expectedFlux[iPoint+fOffset] += p[iPoint];
              expectedFlux[iPoint+fOffset] += scaledPressure[iPoint];
          }
          fOffset += numNodesDomain;
          expectedFlux[iPoint+fOffset] =
            //(rhoE[iPoint]+p[iPoint])*velocity[vectorOffset+iPoint];
            (rhoE[iPoint]+p[iPoint])*velHat[vectorOffset+iPoint];
        }
        int fluxDim = iDim + 1;
        
        FC_MODULE(euler,uniformflux,EULER,UNIFORMFLUX)
          (&numDim,&fluxDim,&numNodes[0],&numNodesDomain,&fortranBufferInterval[0],
           &gridMetric[0],&rho[0],&rhoV[0],&rhoE[0],
           &velHat[0],&p[0],&scaledPressure[0],&eulerFlux[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++){
            if(std::abs(eulerFlux[fOffset+iPoint] - expectedFlux[fOffset+iPoint]) > 1e-15){
              if(!triggered){
                std::cout << "Actual: " << eulerFlux[fOffset+iPoint] << std::endl
                          << "Expected:   " << expectedFlux[fOffset+iPoint] << std::endl;
                triggered = true;
              }
              uniformFlux[fDim] = false;
            }
          }
        }
        for(int fDim = 0;fDim < (numDim+2);fDim++){
          if(!uniformFlux[fDim]){
            dimUniformFlux[numDim-2] = false;
            std::cout << "uniform flux failed in dim,fdim " << iDim 
                      << "," << fDim << std::endl;
          }
        }
        
      }
    }

    plascom2::operators::sbp::base sbp12;
    plascom2::operators::sbp::Initialize(sbp12,2);
    
    // Forticate the stencil starts
    for(int iStencil = 0;iStencil < sbp12.numStencils;iStencil++)
      sbp12.stencilStarts[iStencil]++;
    
    int numStencils = sbp12.numStencils;
    int boundaryDepth = 1;
    
    
    // Create stencil connectivity
    std::vector<int> stencilConnectivity(numDim*numNodesDomain,0);
    plascom2::operators::sbp::CreateStencilConnectivity(numDim,&numNodes[0],&fortranBufferInterval[0],
                                                        boundaryDepth,&stencilConnectivity[0],true);
    
    
    std::vector<size_t> applyNodes(numNodesDomain,0);
    size_t nNode = 1;
    for(size_t iNode = 0;iNode < numNodesDomain;iNode++)
      applyNodes[iNode] = nNode++;
    
    int *stencilSizes = sbp12.stencilSizes;
    int *stencilStarts = sbp12.stencilStarts;
    int numStencilValues = sbp12.numValues;
    double *stencilWeights = sbp12.stencilWeights;
    int *stencilOffsets = sbp12.stencilOffsets;
    int *stencilID = &stencilConnectivity[0];
    int numScalar = 0;
    double *scalarBuffer = NULL;
    
    // Note, RHS *must* be zero upon entry to Euler kernel; Euler kernel
    // only accumulates to RHS vector, doesn't "set" its value.
    std::vector<double> rhsBuffer((numDim+2)*numNodesDomain,0.0);

    FC_MODULE(euler,uniformrhs,EULER,UNIFORMRHS)
      (&numDim,&numNodes[0],&numNodesDomain,&fortranBufferInterval[0],
       &fortranBufferInterval[0],&gridMetric[0],
       &numStencils,&numStencilValues,stencilSizes,stencilStarts,
       stencilOffsets,stencilWeights,stencilID,
       &rho[0],&rhoV[0],&rhoE[0],
       &velHat[0],&p[0],&rhsBuffer[0],&rhsBuffer[numNodesDomain],
       &rhsBuffer[(numDim+1)*numNodesDomain]);
    
    bool eulerRHSWorks = true;
    
    
    for(size_t iNode = 0;iNode < numNodesDomain;iNode++){
      if(std::abs(rhsBuffer[iNode]) > 1e-16){
        std::cout << "drho[" << iNode << "] = " << rhsBuffer[iNode] << std::endl;
        eulerRHSWorks = false;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        if(std::abs(rhsBuffer[iNode+(iDim+1)*numNodesDomain]) > 1e-16){
          std::cout << "dV[" << iDim << "][" << iNode << "] = " 
                    << rhsBuffer[iNode+(iDim+1)*numNodesDomain] 
                    << std::endl;
          eulerRHSWorks = false;
        }
      }
      if(std::abs(rhsBuffer[iNode+(numDim+1)*numNodesDomain]) > 1e-16){
        std::cout << "dE[" << iNode << "] = " 
                  << rhsBuffer[iNode+(numDim+1)*numNodesDomain] << std::endl;
        eulerRHSWorks = false;
      }
    } 
    dimEulerRHS[numDim-2] = eulerRHSWorks;
  }

  bool computeDVWorks = true;
  bool eulerRHSWorks = true;
  bool uniformGridMetricWorks = true;
  bool uniformFluxWorks = true;
  // Make sure these tests worked in 2d and 3d
  for(int nDim = 2;nDim <= 3;nDim++){
    computeDVWorks = computeDVWorks && dimComputeDV[nDim-2];
    eulerRHSWorks  = eulerRHSWorks && dimEulerRHS[nDim-2];
    uniformGridMetricWorks = uniformGridMetricWorks && dimUniformGridMetric[nDim-2];
    uniformFluxWorks = uniformFluxWorks&&dimUniformFlux[nDim-2];
  }
  
  serialUnitResults.UpdateResult("Euler:Kernels:DV",computeDVWorks);
  serialUnitResults.UpdateResult("Euler:Kernels:UniformRHS",eulerRHSWorks);
  serialUnitResults.UpdateResult("Euler:Kernels:UniformFlux",uniformFluxWorks);
  serialUnitResults.UpdateResult("Grid:Kernels:ApplyUniformGridMetric",uniformGridMetricWorks);
};

  
void TestEulerRHS(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{
  std::cout << "TestEulerRHS..." << 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 parTestEulerRHS = true;
  std::vector<size_t> gridSizes(3,4);


  if(testfixtures::CreateSimulationFixtures(simGrid,simHalo,sbpOperator,
                                            gridSizes,2,testComm,&messageStream)){
    parTestEulerRHS = false;
    std::cout << "TestEulerRHS failed in CreateSimulation Fixtures" << std::endl;
    return;
  }

/*
  // Test simple halo exchange to verify Communicators
  std::cout << "Verify Communicators" << std::endl;
  std::vector<int> gridNeighbors;
  simGrid.Communicator().CartNeighbors(gridNeighbors);

  std::vector<double> &GridCoords(simGrid.CoordinateData());
  int messageId = simHalo.CreateMessageBuffers(2); 
  std::cout << "Pack" << std::endl;
  simHalo.PackMessageBuffers(messageId,0,2,&GridCoords[0]);
  std::cout << "Send" << std::endl;
  simHalo.SendMessage(messageId,gridNeighbors,simGrid.Communicator());
  std::cout << "Receive" << std::endl;
  simHalo.ReceiveMessage(messageId,gridNeighbors,simGrid.Communicator());
  std::cout << "Unpack" << std::endl;
  simHalo.UnPackMessageBuffers(messageId,0,2,&GridCoords[0]);
  
  simGrid.Communicator().Barrier();
*/

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

  if(testfixtures::euler::SetupEulerState(simGrid,simState,paramState,0)){
    parTestEulerRHS = false;
    std::cout << "TestEulerRHS failed in SetupEulerState" << std::endl;
    return;
  }

  std::cout << messageStream.str() << std::endl;
  
  size_t numPointsBuffer = simGrid.BufferSize();

  double t = 0.0;
  double gamma = 1.4;
  double dt = .001;
  double cfl = 1.0;
  double Re  = 0.0;
  int flags = 1;

  int numDim = 3;

  paramState.SetFieldBuffer("gamma",&gamma);
  paramState.SetFieldBuffer("inputDT",&dt);
  paramState.SetFieldBuffer("inputCFL",&cfl);
  paramState.SetFieldBuffer("refRe",&Re);
  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,1.0);
  std::vector<double> dRhoV(3*numPointsBuffer,0.0);
  std::vector<double> dRhoE(3*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);


  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 << "Before initialize" << std::endl;

  if(eulerRHS.Initialize(simGrid,simState,paramState,sbpOperator)){
    parTestEulerRHS = false;
    std::cout << "TestEulerRHS failed in eulerRHS.Initialize" << std::endl;
    return;
  }

  std::cout << "Done initialize" << std::endl;
 
  if(eulerRHS.SetupRHS(simGrid,simState,rhsState)){
    parTestEulerRHS = false;
    std::cout << "TestEulerRHS failed in eulerRHS.SetupRHS" << std::endl;
    return;
  }

  std::cout << "Done setupRHS" << std::endl;

  const std::vector<double *> &inputStateBuffers(eulerRHS.StateBuffers());
  const std::vector<double *> &outputStateBuffers(eulerRHS.RHSBuffers());

  
  // check conserved quantities
  if(inputStateBuffers[0] != &rho[0])
    parTestEulerRHS = false;
  for(int iDim = 0;iDim < numDim;iDim++){
    if(inputStateBuffers[iDim+1] != &rhoV[iDim*numPointsBuffer])
      parTestEulerRHS = false;
  }
  if(inputStateBuffers[numDim+1] != &rhoE[0])
    parTestEulerRHS = false;
  
  parTestEulerRHS = CheckResult(parTestEulerRHS,testComm);
  if(!parTestEulerRHS){
    std::cout << "Problem 1" << std::endl;
  }
  // check RHS buffers
  if(outputStateBuffers[0] != &dRho[0])
    parTestEulerRHS = false;
  if(outputStateBuffers[1] != &dRhoV[0])
    parTestEulerRHS = false;
  if(outputStateBuffers[numDim+1] != &dRhoE[0])
    parTestEulerRHS = false;
  parTestEulerRHS = CheckResult(parTestEulerRHS,testComm);
  if(!parTestEulerRHS){
    std::cout << "Problem 2" << std::endl;
  }

  std::cout << "Done Problem 1 and 2" << std::endl;

    // Call SetupRHS again, with switched arguments
  if(eulerRHS.SetupRHS(simGrid,rhsState,simState)){
    parTestEulerRHS = false;
  }
  parTestEulerRHS = CheckResult(parTestEulerRHS,testComm);
  if(!parTestEulerRHS){
    std::cout << "Problem 3" << std::endl;
  }

  // check conserved quantities
  if(inputStateBuffers[0] != &dRho[0])
    parTestEulerRHS = false;
  for(int iDim = 0;iDim < numDim;iDim++)
    if(inputStateBuffers[iDim+1] != &dRhoV[iDim*numPointsBuffer])
      parTestEulerRHS = false;
  if(inputStateBuffers[numDim+1] != &dRhoE[0])
    parTestEulerRHS = false;
  parTestEulerRHS = CheckResult(parTestEulerRHS,testComm);
  if(!parTestEulerRHS){
    std::cout << "Problem 4" << std::endl;
  }

  std::cout << "Done Problem 3 and 4" << std::endl;

  // check RHS buffers
  if(outputStateBuffers[0] != &rho[0])
    parTestEulerRHS = false;
  if(outputStateBuffers[1] != &rhoV[0])
    parTestEulerRHS = false;
  if(outputStateBuffers[numDim+1] != &rhoE[0])
    parTestEulerRHS = false;
  parTestEulerRHS = CheckResult(parTestEulerRHS,testComm);
  if(!parTestEulerRHS){
    std::cout << "Problem 5" << std::endl;
  }

  std::cout << "Done Problem 5" << std::endl;
  
  int returnCode = eulerRHS.RHS(0);
  if(returnCode > 0){
    parTestEulerRHS = false;
  }

  parTestEulerRHS = CheckResult(parTestEulerRHS,testComm);
  if(!parTestEulerRHS){
    std::cout << "Problem 6" << " Return code " << returnCode << std::endl;
  }

  std::cout << "Done Problem 6" << std::endl;
  
  parallelUnitResults.UpdateResult("Euler:RHS:RunsInParallel",parTestEulerRHS);

};

void TestEulerRHS2(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
{
  std::cout << "TestEulerRHS2..." << 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 parTestEulerRHS = true;
  std::vector<size_t> gridSizes(3,4);

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

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

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

  double t = 0.0;
  double gamma = 1.4;
  double dt = .001;
  double cfl = 1.0;
  double Re  = 0.0;
  //  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("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;
  
  if(eulerRHS.Initialize(simGrid,simState,paramState,sbpOperator)){
    parTestEulerRHS = false;
    return;
  }

  if(eulerRHS.SetupRHS(simGrid,simState,rhsState)){
    parTestEulerRHS = false;
    return;
  }

  simGrid.SetNumThreads(1);
  eulerRHS.InitThreadIntervals();
 
  
  int returnCode = eulerRHS.RHS(0);
  if(returnCode > 0){
    parTestEulerRHS = false;
  }

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

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

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

  parallelUnitResults.UpdateResult("Euler:RHS:3dWorksInParallel",parTestEulerRHS);


};