mixedhybridcompositional.cpp

Go to the documentation of this file.

#include "mixedhybridcompositional.h"
#include "numericmethods.h"
#include "hdf5orthowriter.h"
#include "gnuplotanim.h"


MixedHybridCompositional::MixedHybridCompositional(OrthoMesh &mesh,Function3D &fInitP,Function3D &fPrescribedVelocities, Function3D &fPrescribedPression, const VecDouble &K, const VecDouble &vPor,  GeneralFunctionInterface &fMobT,  GeneralFunctionInterface &fFracFlux,double grav, double dt,  FlashCompositional &flash,LinearSolver &solver,int debugLevel,bool bCheckNoBC) 
  :m_mesh(mesh),
   m_fMobT(fMobT),
   m_fFracFlux(fFracFlux),
   m_vPor(vPor),
   m_vK(K),
   m_grav(grav),
   m_dt(dt),
   m_solver(solver),
   m_flash(flash),
   _debugLevel(debugLevel)
{
  
  
  this->setBoundaryCondition(m_bc,m_mesh,fPrescribedVelocities,fPrescribedPression,bCheckNoBC);
  this->buildSparsityPattern(mesh,m_sparsePattern);
  m_A.reinit(m_sparsePattern);
  m_B.reinit(m_mesh.numCells());
  m_sol.reinit(m_mesh.numCells(),NAN);
  m_Vq.reinit(m_mesh.numFaces());
  //Set the initial values for the pressure;
  mesh.setCentralValuesFromFunction(fInitP,m_sol);


  m_vInvH2.reinit(3);
  m_vInvH2(OrthoMesh::NORMAL_X)=1.0/(pow(mesh.get_dx(),2));
  m_vInvH2(OrthoMesh::NORMAL_Y)=1.0/(pow(mesh.get_dy(),2));
  m_vInvH2(OrthoMesh::NORMAL_Z)=1.0/(pow(mesh.get_dz(),2));
  m_inv_hz=1.0/mesh.get_dz();

  
}

 MixedHybridCompositional::~MixedHybridCompositional() 
{
}











void MixedHybridCompositional::getVelocitiesAtFaces(Matrix &vel) 
{
  assert(vel.m() == m_mesh.numFaces() && vel.n() == 3);
  
  vel = 0.0;
  //Get the face iterator and the number of faces.
  OrthoMesh::Face_It face = m_mesh.begin_face();
  OrthoMesh::Face_It endf = m_mesh.end_face();
    
  //For each face
  for (;face!=endf;face++)
  {
    vel(face->index(),face->getNormalOrientation()) = m_Vq(face->index());
  }
}

void MixedHybridCompositional::getVelocitiesAtCells(Matrix &vel) 
{
  m_mesh.projectVelocitiesInFacesToCells(vel,m_Vq);
}



void MixedHybridCompositional::printOutput() 
{
  if (_debugLevel > 0)
  {
    VecDouble v(m_mesh.numCells());
    VecDouble Vv(m_mesh.numVertices());
    HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 

    m_mesh.projectCentralValuesAtVertices(getPressureAtCells(),Vv);
    hdf5.writeScalarField(Vv,"P");



  
    if (_debugLevel > 1)
    {
      Matrix Mvel(m_mesh.numCells(),3);
      FEValuesExtractors::Vector velocities(0);
      this->getVelocitiesAtCells(Mvel);
      if (_debugLevel > 2)
      {
        //printf("\nVelocities at cells' barycenters\n\n");
        //Mvel.print_formatted(std::cout);
      }
      VecDouble v(m_mesh.numCells());
      for (unsigned i=0;i<v.size();i++)
      {
        v(i) = Mvel(i,0); 
      }
      VecDouble minVel(3);
      VecDouble maxVel(3);
      NumericMethods::getMinMaxValueByColumn(Mvel,minVel,maxVel);
      printf("Min Vel <%g,%g,%g> - Max Vel <%g,%g,%g>\n",minVel(0),minVel(1),minVel(2),maxVel(0),maxVel(1),maxVel(2));
  
      m_mesh.projectCentralValuesAtVertices(v,Vv);
      hdf5.writeScalarField(Vv,"Vel");
    }
  }

}




 void MixedHybridCompositional::getNormalVelocityAtFaces(VecDouble &vFNC) 
{
  vFNC=m_Vq;
}


const VecDouble& MixedHybridCompositional::getPressureAtCells()
{
  return m_sol;
}

---