mixedhybridcompositionalfull.cpp

Go to the documentation of this file.

#include "mixedhybridcompositionalfull.h"
#include "numericmethods.h"
#include "hdf5orthowriter.h"


 MixedHybridCompositionalFull::MixedHybridCompositionalFull(OrthoMesh &mesh,Function3D &fPInit,Function3D &fPrescribedVelocities, Function3D &fPrescribedPression, const VecDouble &K, const VecDouble &vPor,Function1D &fMobT,Function1D &fRelMobilities,double grav, double dt,FlashCompositional &flash,LinearSolver &solver,int debugLevel,bool bCheckNoBC) 
   :MixedHybridCompositional(mesh,fPInit,fPrescribedVelocities, fPrescribedPression, K, vPor,fMobT,fRelMobilities,grav,  dt,flash,solver,debugLevel, bCheckNoBC)
{
  m_vDvDm.reinit(mesh.numCells(),2);
  m_vC.reinit(mesh.numCells(),2);
  m_vMobT.reinit(mesh.numCells());
  m_vMobRel.reinit(mesh.numCells());

}




void MixedHybridCompositionalFull::iterate(TransportBase &trans) 
{
  this->assemblySystem(m_A,m_B,m_dt,m_sol,m_vK,m_vPor,m_flash,m_fMobT,m_fFracFlux,m_grav);

  VecDouble PAnt=m_sol;

  m_solver.solve(m_A,m_sol,m_B);

  Point3D p(5,50,50);
  OrthoMesh::Cell_It cell = m_mesh.getCellAt(p);
  //testSolution(cell,m_dt, v1,m_vPor,m_vK,m_sol,PAnt,m_vMobT,v2); 

  //Get the fluxes
  this->getNormalVdtFromPressure(m_Vq,m_mesh,m_bc,m_vK,m_sol,m_vMobT);

}

using NumericMethods::_div;







void MixedHybridCompositionalFull::assemblySystem(SparseMatrix<double> &A,VecDouble &B, double dt, const VecDouble &vP, const VecDouble &vK, const VecDouble &vPor, FlashCompositional &flash, GeneralFunctionInterface &fMobT, GeneralFunctionInterface &fRelMobilities,double grav) 
{
  assert((flash.numPhases() -1) ==  fRelMobilities.getImageDim());
  assert(fMobT.getDomainDim() == flash.numPhases()-1);
  assert(fMobT.getImageDim() == 1);
  assert(fRelMobilities.getDomainDim() == flash.numPhases()-1);

  static VecDouble totalComponents(flash.numComponents());
  static VecDouble dv_dm(flash.numComponents());
  static VecDouble phasesVol(flash.numPhases());
  static VecDouble phasesDensities(flash.numPhases());
  static VecDouble phasesViscosities(flash.numPhases());
  
  double maxStE=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(flash.numPhases()-1);

  
                           

  //ALLOCATE THE flash data
  static FlashData data(flash.numPhases(),flash.numComponents());


  //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
  A=0.0;
  B=0.0;

  //For each cell do
  double BetaG_dt;
  double St;
  for (Index cellIndex=0;cellIndex < nCells; cellIndex++,cell++)
  {
    //Compute/Get the flash data
    flash.flash(cellIndex,data);
    double P=vP(cellIndex);

    flash.getPhasesVolume(P,data,phasesVol);     
    St=phasesVol.l1_norm();


    if (fabs(1-St) > maxStE) maxStE = fabs(1-St); 
    
    BetaG_dt=vPor(cellIndex)*flash.getFluidCompressibility(P,data)/(dt);
    if (BetaG_dt < -0.0)
    {
      printf("Error BetaG_dt = %g\n",BetaG_dt);
      data.print();
      abort();
    }
    A.add(cellIndex,cellIndex,BetaG_dt);
    B(cellIndex)+=BetaG_dt*P+vPor(cellIndex)*(St-1.0)/dt;


                              
    
    
    //Get viscosities
    flash.getPhasesViscosities(P,data,phasesViscosities);
    fMobT.setParameters(phasesViscosities);
    fRelMobilities.setParameters(phasesViscosities);



    
    //Get vector of saturations, remember that the last saturation
    //is a function of the others because the sum of all saturations
    //must be one. So we dont consider the last saturation as a variable
    //for the mobilities
    for (unsigned k=0;k<vS.size();k++)
      vS(k) = phasesVol(k)/St;
    m_vMobT(cellIndex)=fMobT(vS);
    m_vMobRel(cellIndex)=fRelMobilities(vS);

    //Get total conservation components
    data.getTotalComponentsMoles(totalComponents);

    //Get Volume derivatives in terms of the total moles
    flash.getTotalVolumeDerivatives(P,data,dv_dm);

    m_vDvDm(cellIndex,0)=dv_dm(0);
    m_vDvDm(cellIndex,1)=dv_dm(1);
    m_vC(cellIndex,0)=_div(data.getMoles(0,0),phasesVol(0))*fRelMobilities(vS)*fMobT(vS);
    m_vC(cellIndex,1)=( _div(data.getMoles(0,1), phasesVol(0)) * fRelMobilities(vS) +
                        _div(data.getMoles(1,1), phasesVol(1))*(1-fRelMobilities(vS)))*fMobT(vS);
    //data.print();
    //printf("%g %g %g %g\n",phasesVol(0),phasesVol(1),fRelMobilities.value(vS),fMobT.value(vS));
    //printf("%g %g\n", m_vC(cellIndex,0), m_vC(cellIndex,1));
           

  }


  
  //Now with all the information gathered from the flash in each cell
  //build the system

  OrthoMesh::Face_It face = m_mesh.begin_face();
  unsigned nFaces = m_mesh.numFaces();
  



  //For each face
  for (unsigned faceIndex=0;faceIndex < nFaces; faceIndex++,face++)
  {
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);
    
    //Make sure the face is not at boundary
    if (!face->at_boundary())
    {
      double inv_dx2 = m_vInvH2(face->getNormalOrientation());


      double K1 = vK(c1)*m_vC(c1,0);
      double K2 = vK(c2)*m_vC(c2,0);
      double Ceff=2*_div(K1*K2,(K1+K2));
      A.add(c1,c1, +m_vDvDm(c1,0)*Ceff * inv_dx2);
      A.add(c1,c2, -m_vDvDm(c1,0)*Ceff * inv_dx2);
      A.add(c2,c2, +m_vDvDm(c2,0)*Ceff * inv_dx2);
      A.add(c2,c1, -m_vDvDm(c2,0)*Ceff * inv_dx2);


      K1 = vK(c1)*m_vC(c1,1);
      K2 = vK(c2)*m_vC(c2,1);
      Ceff=2*_div(K1*K2,(K1+K2));
      A.add(c1,c1, +m_vDvDm(c1,1)*Ceff * inv_dx2);
      A.add(c1,c2, -m_vDvDm(c1,1)*Ceff * inv_dx2);
      A.add(c2,c2, +m_vDvDm(c2,1)*Ceff * inv_dx2);
      A.add(c2,c1, -m_vDvDm(c2,1)*Ceff * inv_dx2);
    }
  }
  //Now process the boundary conditions
    
  HybridFaceBC::iterator bcEnd = m_bc.end();
  for(HybridFaceBC::iterator bcIt = m_bc.begin();bcIt!=bcEnd;bcIt++)
  {
    OrthoMesh::Face_It face = bcIt->face;
    assert(face->at_boundary());
      
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);

    if (bcIt->bcType == Dirichlet)
    {
      double inv_dx2 =  m_vInvH2(face->getNormalOrientation());
      double c = (c1 != OrthoMesh::invalidIndex()) ? c1 : c2;
      double Coeff = +2*vK(c)*(m_vC(c,0)*m_vDvDm(c,0) + m_vC(c,1)*m_vDvDm(c,1))*inv_dx2;
        
      A.add(c,c, Coeff);
      B(c) += Coeff*bcIt->value;
    }
    else //Neumann condition
    {
      if (bcIt->value != 0.0)
      {
        throw new Exception("The fully compositional model only accept null newmann condition");
                                 
      }
      if (c1 != OrthoMesh::invalidIndex())
        B(c1) += -bcIt->value;
      else
        B(c2) += +bcIt->value;
    }
  }

  HDF5OrthoWriter::getHDF5OrthoWriter().setVariable("St_Error",maxStE);
  maxStE=0;
  
}







void MixedHybridCompositionalFull::getNormalVdtFromPressure(VecDouble &Vq,OrthoMesh &mesh,HybridFaceBC &bc,const VecDouble &K,const VecDouble &vP,const VecDouble &vMobT) 
{
  assert(Vq.size() == mesh.numFaces());
  assert(vP.size() == mesh.numCells());
  assert(vMobT.size() == mesh.numCells());
  assert(K.size() == mesh.numCells());
  
  
  //Get the face iterator and the number of faces.
  OrthoMesh::Face_It face = mesh.begin_face();
  OrthoMesh::Face_It endf = mesh.end_face();

  using NumericMethods::_div;  
  //For each face
  for (;face!=endf;face++)
  {
    if (!face->at_boundary())
    {
      unsigned c1,c2;
      face->getAdjCellIndices(c1,c2);
      double K1 = K(c1)*vMobT(c1);
      double K2 = K(c2)*vMobT(c2);
      double Keff = 2*K1*K2/(K1+K2)*face->areaPerCellVol();
      Vq(face->index()) = - Keff * (vP(c2) - vP(c1));
      //std::cout << vMobT(c1) << " ";
      
      }
  }

  //Now for boundary conditions
    
  HybridFaceBC::iterator bcEnd = bc.end();
  for(HybridFaceBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
  {
    OrthoMesh::Face_It face = bcIt->face;
    assert(face->at_boundary());
      
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);

    if (bcIt->bcType == Dirichlet)
    {
      if (c1 != OrthoMesh::invalidIndex())
      {
        double Keff = 2*K(c1)*vMobT(c1)*face->areaPerCellVol();
        Vq(face->index()) = - Keff*(bcIt->value - vP(c1)); 
        //      printf("Pressions R: %d %g %g = %g\n",c1,bcIt->value, m_sol(c1),Vq(bcIt->faceIndex));
      }
      else
      {
        double Keff = 2*K(c2)*vMobT(c2)*face->areaPerCellVol();
        Vq(face->index()) = Keff*(bcIt->value - vP(c2));
        //printf("Pressions L: %d %g %g = %g\n",c2,bcIt->value, m_sol(c2),Vq(bcIt->faceIndex));
      }
    }
    else
    {
      if (c1 != OrthoMesh::invalidIndex())
        Vq(face->index()) = bcIt->value/face->area();
      else
        Vq(face->index()) = bcIt->value/face->area();
    }
  }

}


 MixedHybridCompositionalFull::~MixedHybridCompositionalFull() 
{

}





void MixedHybridCompositionalFull::getTotalMassFluxFromPressure(VecDouble &Vq,OrthoMesh &mesh,const HybridFaceBC &bc,const VecDouble &K,const VecDouble &vP,const Matrix &MC) 
{
  assert(Vq.size() == mesh.numFaces());
  assert(vP.size() == mesh.numCells());
  assert(K.size() == mesh.numCells());
 
  
  //Get the face iterator and the number of faces.
  OrthoMesh::Face_It face = mesh.begin_face();
  OrthoMesh::Face_It endf = mesh.end_face();

  //  m_vC.print(std::cout,8);
    
  //For each face
  for (;face!=endf;face++)
  {
    if (!face->at_boundary())
    {
      unsigned c1,c2;
      face->getAdjCellIndices(c1,c2);
      //double K1 = K(c1)*MC(c1,0) + K(c1)*MC(c1,1);
      //double K2 = K(c2)*MC(c2,0) + K(c2)*MC(c2,1);
      double K1 = K(c1)*MC(c1,0);
      double K2 = K(c2)*MC(c2,0);
      double K3 = K(c1)*MC(c1,1);
      double K4 = K(c2)*MC(c2,1);
      double Keff = 2*(_div(K1*K2,K1+K2) + _div(K3*K4,K3+K4)) *face->areaPerCellVol();
      Vq(face->index()) = - Keff * (vP(c2) - vP(c1));
    }
  }

  //Now for boundary conditions
    
  HybridFaceBC::const_iterator bcEnd = bc.end();
  for(HybridFaceBC::const_iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++)
  {
    OrthoMesh::Face_It face = bcIt->face;
    assert(face->at_boundary());
      
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);

    if (bcIt->bcType == Dirichlet)
    {
      if (c1 != OrthoMesh::invalidIndex())
      {
        double Keff = 2*K(c1)*(MC(c1,0) + MC(c1,1))*face->areaPerCellVol();
        Vq(face->index()) = - Keff*(bcIt->value - vP(c1)); 
        //      printf("Pressions R: %d %g %g = %g\n",c1,bcIt->value, m_sol(c1),Vq(bcIt->faceIndex));
      }
      else
      {
        double Keff = 2*K(c2)*(MC(c2,0) + MC(c2,1))*face->areaPerCellVol();
        Vq(face->index()) = Keff*(bcIt->value - vP(c2));
        //printf("Pressions L: %d %g %g = %g\n",c2,bcIt->value, m_sol(c2),Vq(bcIt->faceIndex));
      }
    }
    else
    {
      if (c1 != OrthoMesh::invalidIndex())
        Vq(face->index()) = bcIt->value/face->area();
      else
        Vq(face->index()) = bcIt->value/face->area();
    }
  }

  
}





 void MixedHybridCompositionalFull::printOutput() 
{
  MixedHybridCompositional::printOutput();
  
  if (_debugLevel > 1)
  {
    Matrix VelM(m_mesh.numVertices(),3);
    VecDouble Vf(m_mesh.numFaces());
    VecDouble Vv(m_mesh.numVertices());
    VecDouble Vc(m_mesh.numCells());
    HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 
    //getTotalMassFluxFromPressure(Vf,m_mesh,getFaceBC(),m_vK,m_sol,m_vC);
    //m_mesh.projectVelocitiesInFacesToVertices(VelM,Vf);
    //NumericMethods::matrixGetColumn(Vv,VelM,0);
    //hdf5.writeScalarField(Vv,"VelM");

  }
  


  
}


 void MixedHybridCompositionalFull::getNormalMassFluxAtFaces(VecDouble &vFNC) 
{
  getTotalMassFluxFromPressure(vFNC,m_mesh,m_bc,m_vK,m_sol,m_vC);
}

---