mixedhybridincompressible.cpp

Go to the documentation of this file.

#include "mixedhybridincompressible.h"
#include "vecdouble.h"
#include "linearsolver.h"
#include "transportbase.h"
#include "hdf5orthowriter.h"


MixedHybridIncompressible::MixedHybridIncompressible(OrthoMesh &mesh,Function3D &fPrescribedVelocities, Function3D &fPrescribedPression, const VecDouble &K,   GeneralFunctionInterface &fMobT,  GeneralFunctionInterface &fFracFlux,double grav, VecDouble &densities,LinearSolver &solver,int debugLevel) 
  :m_mesh(mesh),
   m_fMobT(fMobT),
   m_fFracFlux(fFracFlux),
   m_densities(densities),
   m_vK(K),
   m_grav(grav),
   m_solver(solver),
   _debugLevel(debugLevel)
{
  
  
  this->setBoundaryCondition(m_bc,m_mesh,fPrescribedVelocities,fPrescribedPression,true);
  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());
  m_vGravTerm.reinit(m_mesh.numCells());
  m_vMobT.reinit(m_mesh.numCells());
  

}


 void MixedHybridIncompressible::iterate(TransportBase &trans) 
{
  assert(m_fFracFlux.getImageDim() == 1);
  assert(m_fMobT.getDomainDim() == 1);
  assert(m_fMobT.getImageDim() == 1);
  assert(m_fFracFlux.getDomainDim() == 1);
  assert(m_grav >= 0);
  
  //Allocate the vector of saturations. Since the sum
  //must be one the saturation of the last phase is
  //a function of the others, so we just have numPhases()-1
  //Independent saturation variables
  static VecDouble vS(1);




  //Get the cell iterator and the number of cells.
  OrthoMesh::Cell_It cell = m_mesh.begin_cell();
  Index nCells = m_mesh.numCells();

  //Write Zeros in the system
  m_A=0.0;
  m_B=0.0;

  ArrayOfVecDouble &transSol = trans.getSolutionAtCells();
  //For each cell do
  for (Index cellIndex=0;cellIndex < nCells; cellIndex++,cell++)
  {
    vS(0)=transSol(cellIndex,0);
    //transSol.getVecValues(cellIndex,&vS);

    m_vMobT(cellIndex)=m_fMobT(vS);
    
    

    double dMob = m_fFracFlux(vS,0);
    double mobAcum=dMob;
    double dGravTerm;
    dGravTerm=dMob*m_densities(0);
    for (unsigned k=1;k<vS.size();k++)
    {
      dMob = m_fFracFlux(vS,k);
      dGravTerm+=dMob*m_densities(k);
      mobAcum+=dMob;  
    }
    dGravTerm+=(1.0-mobAcum)*m_densities(vS.size());

    //Remember that the gravity is in the opposite direction
    //of the Z axis (grav acceleration vector  = <0,0,-grav>)
    //That is why we use the minus sign  here
    m_vGravTerm(cellIndex)=-dGravTerm*m_grav;

    assert(mobAcum<=1.0000000001);
  }
  MixedHybridBase::assemblySystem(m_A,m_B,m_mesh,m_bc,m_vK,m_vMobT,m_vGravTerm);

  //Solve the system, storing the solution on m_sol.
  m_solver.solve(m_A,m_sol,m_B);

  
  //Update the fluxes
  //this->getNormalVelocitiesFromPressure(m_Vq,m_mesh,m_bc,m_vK,m_sol,m_vMobT,m_vGravTerm);

}


 void MixedHybridIncompressible::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 >= 3)
    {
      VecDouble vF(m_mesh.numFaces());
      getNormalVelocityAtFaces(vF);
      Matrix M(m_mesh.numVertices(),3);
      m_mesh.projectVelocitiesInFacesToVertices(M,vF);
      M.copyColumn(Vv,2);
      hdf5.writeScalarField(Vv,"VelZ");
      M.copyColumn(Vv,0);
      hdf5.writeScalarField(Vv,"VelX");
      M.copyColumn(Vv,1);
      hdf5.writeScalarField(Vv,"VelY");

      
    }

  }
  
}

---