TestEOS.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;
typedef simulation::grid::parallel_blockstructured grid_t;

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

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

  std::vector<bool> dimComputeEOS(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;
    
    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;
    }
    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;
    double v = vVec[0];

    std::vector<double> expectedPressure(numNodesDomain,1.0);
    std::vector<double> expectedTemperature(numNodesDomain,1.0);
    
    std::vector<double> ke(numNodesDomain,0.0);
    std::vector<double> rho(numNodesDomain,1.0);
    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];
      }
    } 
    // 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.);
    std::vector<double> p(numNodesDomain,0.);
    std::vector<double> rhom1(numNodesDomain,0.);
    std::vector<double> velocity(numDim*numNodesDomain,0.);
    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++)
      myDVBuffers[iDim+3] = &velocity[iDim*numNodesDomain];
    myDVBuffers[numDim+3] = &internalEnergy[0];

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

    //std::cout << "bufferExtent " << bufferExtent << std::endl;
    //std::cout << "bufferInterval " << bufferInterval << std::endl;

    // 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;
    }

    //std::cout << "Before ComputeDVBuffer" << std::endl;
    
    // first calculate the specific volume, velocity, and internal energy
    // this is testing elsewhere
    if(euler::util::ComputeDVBuffer(bufferExtent,numNodes,myStateBuffers,myDVBuffers)) {
      std::cout << "ComputeDVBuffer failed." << std::endl;
    }

    //std::cout << "Before ComputePressureBuffer" << 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 computeEOSWorks = true;
    bool pfailed = false;
    bool tfailed = false;

    bool printPressure=true;
    bool printTemperature=true;
    for(int iNode = 0;iNode < numNodesDomain;iNode++){
      if(std::abs(p[iNode] - expectedPressure[iNode])    > 1e-15){
        if(printPressure){
          std::cout << "In EOS_ComputeBuffers, Expected p=" << expectedPressure[iNode]
                    << " got p=" << p[iNode] << " at iNode=" << iNode << std::endl;
          std::cout << "No further messages will be generated for " << std::endl;
          printPressure = false;
        }
        computeEOSWorks = false;
        pfailed = true;
      }

      if(std::abs(T[iNode] - expectedTemperature[iNode]) > 1e-15){
        if(printTemperature){
          std::cout << "In EOS_ComputeBuffers, Expected T=" << expectedTemperature[iNode]
                    << " got T=" << T[iNode] << " at iNode=" << iNode << std::endl;
          std::cout << "No further messages will be generated for " << std::endl;
          printTemperature = false;
        }
        computeEOSWorks = false;
        tfailed = true;
      }
    }
    
    dimComputeEOS[numDim-2] = computeEOSWorks ;
    if(pfailed){
      std::cout << "Pressure failed" << std::endl;
    }
    if(tfailed){
      std::cout << "Temperature failed" << std::endl;
    }
  }

  bool computeEOSWorks = true;
  // Make sure these tests worked in 2d and 3d
  for(int nDim = 2;nDim <= 3;nDim++){
    computeEOSWorks = computeEOSWorks && dimComputeEOS[nDim-2];
  }
  
  serialUnitResults.UpdateResult("EOS:ComputeBuffers",computeEOSWorks);
};