mixedhybridcompositionalsimplepw.cpp

Go to the documentation of this file.

#include "mixedhybridcompositionalsimplepw.h"
#include "hdf5orthowriter.h"
#include "numericmethods.h"


MixedHybridCompositionalSimplePw::MixedHybridCompositionalSimplePw(OrthoMesh &mesh,Function3D &fInitP,Function3D &fPrescribedVelocities, Function3D &fPrescribedPressure, Function3D &fTransportBC,const VecDouble &K, const VecDouble &vPor,  GeneralFunctionInterface &fMobT,  GeneralFunctionInterface &fFracFlux,GeneralFunctionInterface &pc,double grav, double dt,  FlashCompositional &flash,LinearSolver &solver,int debugLevel)
 :MixedHybridCompositional(mesh,fInitP,fPrescribedVelocities, fPrescribedPressure, K, vPor,fMobT,fFracFlux,grav,  dt,flash,solver,debugLevel, false),m_fpc(pc)
{
  m_vGravTerm.reinit(m_mesh.numCells());
  m_vMobT.reinit(m_mesh.numCells());
  m_vPc.reinit(m_mesh.numCells());
  m_vKPc.reinit(m_mesh.numCells());
  m_bFirstIteration=true;
  m_solAnt=m_sol;
  setScalarBC(mesh,m_pcbc,fPrescribedPressure,fTransportBC,pc,flash);
  
}





void MixedHybridCompositionalSimplePw::assemblySystem(SparseMatrix<double> &A,VecDouble &B, double dt, const VecDouble &vPw, const VecDouble &vK, const VecDouble &vPor, FlashCompositional &flash, GeneralFunctionInterface &fMobT, GeneralFunctionInterface &fFracFlux,GeneralFunctionInterface &fPc, double grav) 
{
  assert(flash.numPhases()==2);
  assert((flash.numPhases() -1) ==  fFracFlux.getImageDim());
  assert(fMobT.getDomainDim() == flash.numPhases()-1);
  assert(fMobT.getImageDim() == 1);
  assert(fFracFlux.getDomainDim() == flash.numPhases()-1);
  assert(grav >= 0);
  
  static VecDouble phasesVol(flash.numPhases());
  static VecDouble phasesDensities(flash.numPhases());
  static VecDouble phasesViscosities(flash.numPhases());
  double maxStE=0.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();
  double cellVol=cell->volume();  

  //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=vPw(cellIndex);

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


    if (fabs(1-St) > maxStE) maxStE = fabs(1-St); 
    
    BetaG_dt=vPor(cellIndex)*m_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*cellVol);


    B(cellIndex)+=(BetaG_dt*vPw(cellIndex)+vPor(cellIndex)*(St-1.0)/dt)*cellVol;
    

    //B(cellIndex)+=(BetaG_dt*vP(cellIndex))*cellVol;

    //Get viscosities
    //m_flash.getPhasesViscosities(P,data,phasesViscosities);
    //fMobT.setParameters(phasesViscosities);
    //fFracFlux.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;
    double mobT = fMobT(vS);
    double dMob = fFracFlux(vS,0);

    //Total mobility
    m_vMobT(cellIndex)=mobT;


    //Get cappilar pressure value.
    m_vPc(cellIndex)=m_fpc(vS);
    
    //Get the hydraulic permeability of the pc term
    //Compute  K (1-mob(Sw)) = K ft * fo 
    m_vKPc(cellIndex)=mobT*(1-dMob)*m_vK(cellIndex);
    
    
    //calculate the phases densities
    const VecDouble& compMass=flash.getComponentsMolarMass();//Get components molar masses
    data.getPhasesDensities(phasesDensities,compMass,phasesVol);

    //Calculate the fractional flow functions
    double mobAcum=dMob;
    double dGravTerm;
    dGravTerm=dMob*phasesDensities(0);
    for (unsigned k=1;k<vS.size();k++)
    {
      dMob = fFracFlux(vS,k);
      dGravTerm+=dMob*phasesDensities(k);
      mobAcum+=dMob;  
    }
    dGravTerm+=(1.0-mobAcum)*phasesDensities(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*grav;


    
    assert(mobAcum<=1.0);

    
  }
  MixedHybridBase::assemblySystem(A,B,m_mesh,m_bc,vK,m_vMobT,m_vGravTerm);
  MixedHybridBase::addPcTerm(B,m_mesh,m_pcbc,m_vKPc,m_vPc);

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




  
}









void MixedHybridCompositionalSimplePw::setScalarBC(OrthoMesh &mesh,ScalarBC &bc,Function3D &fPressure,Function3D &fTransBC,GeneralFunctionInterface &fPc,FlashCompositional &flash) 
{
  assert(flash.numPhases() == 2);
  assert(fTransBC.getImageDim() == flash.numComponents());
  unsigned ncmps = flash.numComponents();
  OrthoMesh::Face_It face = mesh.begin_face();
  OrthoMesh::Face_It endf = mesh.end_face();
  
  VecDouble cmps(flash.numComponents());
  VecDouble phasesVol(flash.numPhases());
  _ScalarBC pcbc;
  VecDouble vS(1);
  FlashData flashData(flash.numPhases(),flash.numComponents());
  for (;face!=endf;face++)
  {
    if (face->at_boundary())
    {
      Point3D p = face->barycenter();
      if (fTransBC.isInDomain(p))
      {
        if (fPressure.isInDomain(p))
        {
          for (unsigned i=0;i<ncmps;i++)
          {
            cmps(i)=fTransBC(p,i);
          }
          double P = fPressure(p);
          flash.flash(P,cmps,flashData);
          flash.getPhasesVolume(P,flashData,phasesVol);
          double factor=1.0/NumericMethods::vectorSum(phasesVol);
          cmps*=factor;
          vS(0)=phasesVol(0)/(phasesVol(0) + phasesVol(1));
          pcbc.value=fPc(vS);
          pcbc.face=face;
          bc.push_back(pcbc);
        }
      }
    }
  }
}


void MixedHybridCompositionalSimplePw::iterate(TransportBase &trans) 
{
  m_solAnt=m_sol;
  this->assemblySystem(m_A,m_B,m_dt,m_sol,m_vK,m_vPor,m_flash,m_fMobT,m_fFracFlux,m_fpc,m_grav);
  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->getNormalVelocitiesFromPressure(m_Vq,m_mesh,m_bc,m_vK,m_sol,m_vMobT,m_vGravTerm,m_pcbc,m_vKPc,m_vPc);

}



void MixedHybridCompositionalSimplePw::rollback()
{
  m_sol=m_solAnt;
}

---