mixedhybridcompositionalsimple.cpp

Go to the documentation of this file.

#include "mixedhybridcompositionalsimple.h"
#include "numericmethods.h"
#include "hdf5orthowriter.h"
#include "hdf5orthowriter.h"
#include "solvergpu_agm.h"
#include <fstream>

 MixedHybridCompositionalSimple::MixedHybridCompositionalSimple(OrthoMesh &mesh,Function3D &fInitP,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,fInitP,fPrescribedVelocities, fPrescribedPression, K, vPor,fMobT,fRelMobilities,grav,  dt,flash,solver,debugLevel, bCheckNoBC)
      {
  m_vGravTerm.reinit(m_mesh.numCells());
  m_vMobT.reinit(m_mesh.numCells());
  v1.reinit(m_mesh.numCells());
  v2.reinit(m_mesh.numCells());
  m_solAnt = m_sol;

 }






void MixedHybridCompositionalSimple::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_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);
  
}





double MixedHybridCompositionalSimple::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()==2);
  assert((flash.numPhases() -1) ==  fRelMobilities.getImageDim());
  assert(fMobT.getDomainDim() == flash.numPhases()-1);
  assert(fMobT.getImageDim() == 1);
  assert(fRelMobilities.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=vP(cellIndex);

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

    v2(cellIndex)=St;

    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);
    v1(cellIndex)=-BetaG_dt*cellVol;


    B(cellIndex)+=(BetaG_dt*vP(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);
    //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);
    
    
    //calculate the phases densities
    const VecDouble& compMass=flash.getComponentsMolarMass();//Get components molar masses
    data.getPhasesDensities(phasesDensities,compMass,phasesVol);

    double dMob = fRelMobilities(vS,0);
    double mobAcum=dMob;
    double dGravTerm;
    dGravTerm=dMob*phasesDensities(0);
    for (unsigned k=1;k<vS.size();k++)
    {
      dMob = fRelMobilities(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);

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


void MixedHybridCompositionalSimple::testSolution(OrthoMesh::Cell_It &cell, double dt, VecDouble& Bg,const VecDouble &por,const VecDouble &K,const VecDouble &P,const VecDouble &vPprev,const VecDouble &vMobT,const VecDouble &St) 
{
  double U = cell->index_up();
  double D = cell->index_down();
  double L = cell->index_left();
  double R  = cell->index_right();
  double F = cell->index_front();
  double B = cell->index_back();
  double I = cell->index();

  double Kr = 2*vMobT(R)*K(R)*vMobT(I)*K(I)/(vMobT(R)*K(R) + vMobT(I)*K(I));
  double Kl = 2*vMobT(L)*K(L)*vMobT(I)*K(I)/(vMobT(L)*K(L) + vMobT(I)*K(I));
  double Ku = 2*vMobT(U)*K(U)*vMobT(I)*K(I)/(vMobT(U)*K(U) + vMobT(I)*K(I));
  double Kd = 2*vMobT(D)*K(D)*vMobT(I)*K(I)/(vMobT(D)*K(D) + vMobT(I)*K(I));
  double Kf = 2*vMobT(F)*K(F)*vMobT(I)*K(I)/(vMobT(F)*K(F) + vMobT(I)*K(I));
  double Kb = 2*vMobT(B)*K(B)*vMobT(I)*K(I)/(vMobT(B)*K(B) + vMobT(I)*K(I));
  Point3D DX = m_mesh.getDX();
  double Pant = vPprev(I);
  double result;
  result=  por(I)*Bg(I)*P(I)/dt +
    (Kr*(P(I)-P(R)) + Kl*(P(I)-P(L)))/(DX(0)*DX(0)) +
    (Ku*(P(I)-P(U)) + Kd*(P(I)-P(D)))/(DX(1)*DX(1)) +
    (Kf*(P(I)-P(F)) + Kb*(P(I)-P(B)))/(DX(2)*DX(2))
    -
    (por(I)*Bg(I)*Pant/dt + por(I)*(St(I)-1.0)/dt);
    
  printf("Residual = %g\n", result);
    

  
}


void MixedHybridCompositionalSimple::assemblyInitialSystem(SparseMatrix<double> &A,VecDouble &B,  HybridFaceBC &bc,const VecDouble &vP, const VecDouble &vK,  FlashCompositional &flash, GeneralFunctionInterface &fMobT, GeneralFunctionInterface &fRelMobilities,double grav) 
{
  assert(flash.numPhases()==2);
  assert((flash.numPhases() -1) ==  fRelMobilities.getImageDim());
  assert(fMobT.getDomainDim() == flash.numPhases()-1);
  assert(fMobT.getImageDim() == 1);
  assert(fRelMobilities.getDomainDim() == flash.numPhases()-1);
  assert(grav >= 0);
  
  VecDouble phasesVol(flash.numPhases());
  VecDouble phasesDensities(flash.numPhases());
  VecDouble phasesViscosities(flash.numPhases());
  
  //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
  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 (Index cellIndex=0;cellIndex < nCells; cellIndex++,cell++)
  {
    //Compute/Get the flash data
    flash.flash(cellIndex,data);
    double P=vP(cellIndex);

    

    //Get viscosities
    //m_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
    flash.getPhasesVolume(P,data,phasesVol);     
    double St=NumericMethods::vectorSum(phasesVol);
    for (unsigned k=0;k<vS.size();k++)
      vS(k) = phasesVol(k)/St;
    m_vMobT(cellIndex)=fMobT(vS);
    
    
    //calculate the phases densities
    const VecDouble& compMass=flash.getComponentsMolarMass();//Get components molar masses
    data.getPhasesDensities(phasesDensities,compMass,phasesVol);

    double dMob = fRelMobilities(vS,0);
    double mobAcum=dMob;
    double dGravTerm;
    dGravTerm=dMob*phasesDensities(0);
    for (unsigned k=1;k<vS.size();k++)
    {
      dMob = fRelMobilities(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,bc,vK,m_vMobT,m_vGravTerm);
}


void MixedHybridCompositionalSimple::findInitialPressure(double referencePressure,double tol) 
{
  
  HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter();
  VecDouble PAnt = m_sol;
  VecDouble vv(m_mesh.numVertices());
  m_flash.execute();
  this->assemblyInitialSystem(m_A,m_B,m_bc,m_sol,m_vK,m_flash,m_fMobT,m_fFracFlux,0.0);
  m_solver.solve(m_A,m_sol,m_B);
  m_flash.execute();
  this->assemblyInitialSystem(m_A,m_B,m_bc,m_sol,m_vK,m_flash,m_fMobT,m_fFracFlux,m_grav);
  m_solver.solve(m_A,m_sol,m_B);
  double error;
  m_mesh.projectCentralValuesAtVertices(m_sol,vv);
  hdf5.writeScalarField(vv,"P0");
  printf("Find Initial Pressure");
  while((error=NumericMethods::vectorMaxRelDiff(PAnt,m_sol)) > tol)
  {
    printf("Error %g\n",error);
    m_flash.execute();
    this->assemblyInitialSystem(m_A,m_B,m_bc,m_sol,m_vK,m_flash,m_fMobT,m_fFracFlux,m_grav);
    PAnt=m_sol;
    m_solver.solveAgain(m_A,m_sol,m_B);
    hdf5.writeScalarField(vv,"P0");
  }
  printf("Error %g\n",error);

}


 void MixedHybridCompositionalSimple::printOutput() 
{
  MixedHybridCompositional::printOutput();

 HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 
  
  VecDouble Vv(m_mesh.numVertices());


  if (_debugLevel >= 3)
  {
    Matrix M(m_mesh.numVertices(),3);
    m_mesh.projectVelocitiesInFacesToVertices(M,m_Vq);
    M.copyColumn(Vv,2);
    hdf5.writeScalarField(Vv,"VelZ");
  }

  m_mesh.projectCentralValuesAtVertices(v2,Vv);
  hdf5.writeScalarField(Vv,"St");

}


 void MixedHybridCompositionalSimple::rollback() 
{
  printf("ROLLLLLLLLLLLLLLLLLL\n");
  m_sol=m_solAnt;
}

---