MixedHybridCompositionalSimple Class Reference

#include <mixedhybridcompositionalsimple.h>

Inheritance diagram for MixedHybridCompositionalSimple:

Collaboration diagram for MixedHybridCompositionalSimple:

Public Member Functions
	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=true)
void	findInitialPressure (double referencePressure=100e+5, double tol=1e-5)
void	testSolution (OrthoMesh::Cell_It &cell, double dt, VecDouble &Bt, const VecDouble &por, const VecDouble &K, const VecDouble &P, const VecDouble &vPprev, const VecDouble &vMobT, const VecDouble &St)
	~MixedHybridCompositionalSimple ()
virtual void	iterate (TransportBase &trans)
virtual void	printOutput ()
virtual void	rollback ()
Protected Member Functions
double	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)
void	assemblyInitialSystem (SparseMatrix< double > &A, VecDouble &B, HybridFaceBC &bc, const VecDouble &vP, const VecDouble &vK, FlashCompositional &flash, GeneralFunctionInterface &fMobT, GeneralFunctionInterface &fRelMobilities, double grav)
Protected Attributes
VecDouble	v1
VecDouble	v2
VecDouble	m_vMobT
VecDouble	m_vGravTerm
Private Attributes
double	_maxStError

Detailed Description

This is the MixedHybridCompositional for the Simple model where we use the Kuhn Tucker Conditions and some other assumptions to describe the compositional model as

B * dp/dt + div(vdt) = por (St -1)/dt

Where B is the compressibility, St is the total saturation and vdt is the total darcy velocity defined in terms of the pressure as

vdt = - mobT K (Grad(P) - mobGravT * grav* <0,0,-1>);

Where grav is the absolute value of the gravity. Note that we assume that the gravity vector is pointing always down in the opposite direction of the the Z axis.

mobGravT is given by Sum (rho_i mob_i) where i=1....number of phases.

Definition at line 27 of file mixedhybridcompositionalsimple.h.

Constructor & Destructor Documentation

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 = `true`
	)

Definition at line 8 of file mixedhybridcompositionalsimple.cpp.

00009     :MixedHybridCompositional(mesh,fInitP,fPrescribedVelocities, fPrescribedPression, K, vPor,fMobT,fRelMobilities,grav,  dt,flash,solver,debugLevel, bCheckNoBC)
00010       {
00011   m_vGravTerm.reinit(m_mesh.numCells());
00012   m_vMobT.reinit(m_mesh.numCells());
00013   v1.reinit(m_mesh.numCells());
00014   v2.reinit(m_mesh.numCells());
00015   m_solAnt = m_sol;
00016 
00017  }

MixedHybridCompositionalSimple::~MixedHybridCompositionalSimple ( ) [inline]

Definition at line 53 of file mixedhybridcompositionalsimple.h.

00053 {}

Member Function Documentation

void MixedHybridCompositionalSimple::assemblyInitialSystem	(	SparseMatrix< double > &	A,
		VecDouble &	B,
		HybridFaceBC &	bc,
		const VecDouble &	vP,
		const VecDouble &	vK,
		FlashCompositional &	flash,
		GeneralFunctionInterface &	fMobT,
		GeneralFunctionInterface &	fRelMobilities,
		double	grav
	)			`[protected]`

Definition at line 184 of file mixedhybridcompositionalsimple.cpp.

00185 {
00186   assert(flash.numPhases()==2);
00187   assert((flash.numPhases() -1) ==  fRelMobilities.getImageDim());
00188   assert(fMobT.getDomainDim() == flash.numPhases()-1);
00189   assert(fMobT.getImageDim() == 1);
00190   assert(fRelMobilities.getDomainDim() == flash.numPhases()-1);
00191   assert(grav >= 0);
00192   
00193   VecDouble phasesVol(flash.numPhases());
00194   VecDouble phasesDensities(flash.numPhases());
00195   VecDouble phasesViscosities(flash.numPhases());
00196   
00197   //Allocate the vector of saturations. Since the sum
00198   //must be one the saturation of the last phase is
00199   //a function of the others, so we just have numPhases()-1
00200   //Independent saturation variables
00201   static VecDouble vS(flash.numPhases()-1);
00202 
00203 
00204   //Allocate the flash data
00205   FlashData data(flash.numPhases(),flash.numComponents());
00206 
00207 
00208   //Get the cell iterator and the number of cells.
00209   OrthoMesh::Cell_It cell = m_mesh.begin_cell();
00210   Index nCells = m_mesh.numCells();
00211 
00212   //Write Zeros in the system
00213   A=0.0;
00214   B=0.0;
00215 
00216   for (Index cellIndex=0;cellIndex < nCells; cellIndex++,cell++)
00217   {
00218     //Compute/Get the flash data
00219     flash.flash(cellIndex,data);
00220     double P=vP(cellIndex);
00221 
00222     
00223 
00224     //Get viscosities
00225     //m_flash.getPhasesViscosities(P,data,phasesViscosities);
00226     //fMobT.setParameters(phasesViscosities);
00227     //fRelMobilities.setParameters(phasesViscosities);
00228 
00229     //Get vector of saturations, remember that the last saturation
00230     //is a function of the others because the sum of all saturations
00231     //must be one. So we dont consider the last saturation as a variable
00232     //for the mobilities
00233     flash.getPhasesVolume(P,data,phasesVol);     
00234     double St=NumericMethods::vectorSum(phasesVol);
00235     for (unsigned k=0;k<vS.size();k++)
00236       vS(k) = phasesVol(k)/St;
00237     m_vMobT(cellIndex)=fMobT(vS);
00238     
00239     
00240     //calculate the phases densities
00241     const VecDouble& compMass=flash.getComponentsMolarMass();//Get components molar masses
00242     data.getPhasesDensities(phasesDensities,compMass,phasesVol);
00243 
00244     double dMob = fRelMobilities(vS,0);
00245     double mobAcum=dMob;
00246     double dGravTerm;
00247     dGravTerm=dMob*phasesDensities(0);
00248     for (unsigned k=1;k<vS.size();k++)
00249     {
00250       dMob = fRelMobilities(vS,k);
00251       dGravTerm+=dMob*phasesDensities(k);
00252       mobAcum+=dMob;  
00253     }
00254     dGravTerm+=(1.0-mobAcum)*phasesDensities(vS.size());
00255 
00256     //Remember that the gravity is in the opposite direction
00257     //of the Z axis (grav acceleration vector  = <0,0,-grav>)
00258     //That is why we use the minus sign  here
00259     m_vGravTerm(cellIndex)=-dGravTerm*grav;
00260 
00261     assert(mobAcum<=1.0);
00262   }
00263 
00264   
00265   MixedHybridBase::assemblySystem(A,B,m_mesh,bc,vK,m_vMobT,m_vGravTerm);
00266 }

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
	)			`[protected]`

Definition at line 46 of file mixedhybridcompositionalsimple.cpp.

00047 {
00048   assert(flash.numPhases()==2);
00049   assert((flash.numPhases() -1) ==  fRelMobilities.getImageDim());
00050   assert(fMobT.getDomainDim() == flash.numPhases()-1);
00051   assert(fMobT.getImageDim() == 1);
00052   assert(fRelMobilities.getDomainDim() == flash.numPhases()-1);
00053   assert(grav >= 0);
00054   
00055   static VecDouble phasesVol(flash.numPhases());
00056   static VecDouble phasesDensities(flash.numPhases());
00057   static VecDouble phasesViscosities(flash.numPhases());
00058   double maxStE=0.0;
00059   //Allocate the vector of saturations. Since the sum
00060   //must be one the saturation of the last phase is
00061   //a function of the others, so we just have numPhases()-1
00062   //Independent saturation variables
00063   static VecDouble vS(flash.numPhases()-1);
00064 
00065 
00066   //Allocate the flash data
00067   static FlashData data(flash.numPhases(),flash.numComponents());
00068 
00069 
00070   //Get the cell iterator and the number of cells.
00071   OrthoMesh::Cell_It cell = m_mesh.begin_cell();
00072   Index nCells = m_mesh.numCells();
00073   double cellVol=cell->volume();  
00074   //Write Zeros in the system
00075   A=0.0;
00076   B=0.0;
00077 
00078   //For each cell do
00079   double BetaG_dt;
00080   double St;
00081   for (Index cellIndex=0;cellIndex < nCells; cellIndex++,cell++)
00082   {
00083     //Compute/Get the flash data
00084     flash.flash(cellIndex,data);
00085     double P=vP(cellIndex);
00086 
00087     flash.getPhasesVolume(P,data,phasesVol);     
00088     St=phasesVol.l1_norm();
00089 
00090     v2(cellIndex)=St;
00091 
00092     if (fabs(1-St) > maxStE) maxStE = fabs(1-St); 
00093     
00094     BetaG_dt=vPor(cellIndex)*m_flash.getFluidCompressibility(P,data)/(dt);
00095     if (BetaG_dt < -0.0)
00096     {
00097       printf("Error BetaG_dt = %g\n",BetaG_dt);
00098       data.print();
00099       abort();
00100     }
00101     A.add(cellIndex,cellIndex,BetaG_dt*cellVol);
00102     v1(cellIndex)=-BetaG_dt*cellVol;
00103 
00104 
00105     B(cellIndex)+=(BetaG_dt*vP(cellIndex)+vPor(cellIndex)*(St-1.0)/dt)*cellVol;
00106     //B(cellIndex)+=(BetaG_dt*vP(cellIndex))*cellVol;
00107 
00108     //Get viscosities
00109     //m_flash.getPhasesViscosities(P,data,phasesViscosities);
00110     //fMobT.setParameters(phasesViscosities);
00111     //fRelMobilities.setParameters(phasesViscosities);
00112     //Get vector of saturations, remember that the last saturation
00113     //is a function of the others because the sum of all saturations
00114     //must be one. So we dont consider the last saturation as a variable
00115     //for the mobilities
00116     for (unsigned k=0;k<vS.size();k++)
00117       vS(k) = phasesVol(k)/St;
00118     m_vMobT(cellIndex)=fMobT(vS);
00119     
00120     
00121     //calculate the phases densities
00122     const VecDouble& compMass=flash.getComponentsMolarMass();//Get components molar masses
00123     data.getPhasesDensities(phasesDensities,compMass,phasesVol);
00124 
00125     double dMob = fRelMobilities(vS,0);
00126     double mobAcum=dMob;
00127     double dGravTerm;
00128     dGravTerm=dMob*phasesDensities(0);
00129     for (unsigned k=1;k<vS.size();k++)
00130     {
00131       dMob = fRelMobilities(vS,k);
00132       dGravTerm+=dMob*phasesDensities(k);
00133       mobAcum+=dMob;  
00134     }
00135     dGravTerm+=(1.0-mobAcum)*phasesDensities(vS.size());
00136 
00137     //Remember that the gravity is in the opposite direction
00138     //of the Z axis (grav acceleration vector  = <0,0,-grav>)
00139     //That is why we use the minus sign  here
00140     m_vGravTerm(cellIndex)=-dGravTerm*grav;
00141 
00142     assert(mobAcum<=1.0);
00143   }
00144   MixedHybridBase::assemblySystem(A,B,m_mesh,m_bc,vK,m_vMobT,m_vGravTerm);
00145 
00146   HDF5OrthoWriter::getHDF5OrthoWriter().setVariable("St_Error",maxStE);
00147   return maxStE;
00148 }

void MixedHybridCompositionalSimple::findInitialPressure	(	double	referencePressure = `100e+5`,
		double	tol = `1e-5`
	)

Definition at line 269 of file mixedhybridcompositionalsimple.cpp.

00270 {
00271   
00272   HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter();
00273   VecDouble PAnt = m_sol;
00274   VecDouble vv(m_mesh.numVertices());
00275   m_flash.execute();
00276   this->assemblyInitialSystem(m_A,m_B,m_bc,m_sol,m_vK,m_flash,m_fMobT,m_fFracFlux,0.0);
00277   m_solver.solve(m_A,m_sol,m_B);
00278   m_flash.execute();
00279   this->assemblyInitialSystem(m_A,m_B,m_bc,m_sol,m_vK,m_flash,m_fMobT,m_fFracFlux,m_grav);
00280   m_solver.solve(m_A,m_sol,m_B);
00281   double error;
00282   m_mesh.projectCentralValuesAtVertices(m_sol,vv);
00283   hdf5.writeScalarField(vv,"P0");
00284   printf("Find Initial Pressure");
00285   while((error=NumericMethods::vectorMaxRelDiff(PAnt,m_sol)) > tol)
00286   {
00287     printf("Error %g\n",error);
00288     m_flash.execute();
00289     this->assemblyInitialSystem(m_A,m_B,m_bc,m_sol,m_vK,m_flash,m_fMobT,m_fFracFlux,m_grav);
00290     PAnt=m_sol;
00291     m_solver.solveAgain(m_A,m_sol,m_B);
00292     hdf5.writeScalarField(vv,"P0");
00293   }
00294   printf("Error %g\n",error);
00295 
00296 }

void MixedHybridCompositionalSimple::iterate ( TransportBase & trans ) [virtual]

Implements DynamicBase.

Reimplemented in MixedHybridCompositionalSimpleFP.

Definition at line 24 of file mixedhybridcompositionalsimple.cpp.

00025 {
00026   m_solAnt=m_sol;
00027   this->assemblySystem(m_A,m_B,m_dt,m_sol,m_vK,m_vPor,m_flash,m_fMobT,m_fFracFlux,m_grav);
00028 
00029   m_solver.solve(m_A,m_sol,m_B);
00030 
00031 
00032   
00033   //Point3D p(5,50,50);
00034   //OrthoMesh::Cell_It cell = m_mesh.getCellAt(p);
00035   //testSolution(cell,m_dt, v1,m_vPor,m_vK,m_sol,PAnt,m_vMobT,v2); 
00036 
00037   //Get the fluxes
00038   this->getNormalVelocitiesFromPressure(m_Vq,m_mesh,m_bc,m_vK,m_sol,m_vMobT,m_vGravTerm);
00039   
00040 }

void MixedHybridCompositionalSimple::printOutput ( ) [virtual]

Reimplemented from MixedHybridCompositional.

Definition at line 299 of file mixedhybridcompositionalsimple.cpp.

00300 {
00301   MixedHybridCompositional::printOutput();
00302 
00303  HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 
00304   
00305   VecDouble Vv(m_mesh.numVertices());
00306 
00307 
00308   if (_debugLevel >= 3)
00309   {
00310     Matrix M(m_mesh.numVertices(),3);
00311     m_mesh.projectVelocitiesInFacesToVertices(M,m_Vq);
00312     M.copyColumn(Vv,2);
00313     hdf5.writeScalarField(Vv,"VelZ");
00314   }
00315 
00316   m_mesh.projectCentralValuesAtVertices(v2,Vv);
00317   hdf5.writeScalarField(Vv,"St");
00318 
00319 }

void MixedHybridCompositionalSimple::rollback ( ) [virtual]

Reimplemented from CompressibleDynamic.

Definition at line 322 of file mixedhybridcompositionalsimple.cpp.

00323 {
00324   printf("ROLLLLLLLLLLLLLLLLLL\n");
00325   m_sol=m_solAnt;
00326 }

void MixedHybridCompositionalSimple::testSolution	(	OrthoMesh::Cell_It &	cell,
		double	dt,
		VecDouble &	Bt,
		const VecDouble &	por,
		const VecDouble &	K,
		const VecDouble &	P,
		const VecDouble &	vPprev,
		const VecDouble &	vMobT,
		const VecDouble &	St
	)

Definition at line 151 of file mixedhybridcompositionalsimple.cpp.

00152 {
00153   double U = cell->index_up();
00154   double D = cell->index_down();
00155   double L = cell->index_left();
00156   double R  = cell->index_right();
00157   double F = cell->index_front();
00158   double B = cell->index_back();
00159   double I = cell->index();
00160 
00161   double Kr = 2*vMobT(R)*K(R)*vMobT(I)*K(I)/(vMobT(R)*K(R) + vMobT(I)*K(I));
00162   double Kl = 2*vMobT(L)*K(L)*vMobT(I)*K(I)/(vMobT(L)*K(L) + vMobT(I)*K(I));
00163   double Ku = 2*vMobT(U)*K(U)*vMobT(I)*K(I)/(vMobT(U)*K(U) + vMobT(I)*K(I));
00164   double Kd = 2*vMobT(D)*K(D)*vMobT(I)*K(I)/(vMobT(D)*K(D) + vMobT(I)*K(I));
00165   double Kf = 2*vMobT(F)*K(F)*vMobT(I)*K(I)/(vMobT(F)*K(F) + vMobT(I)*K(I));
00166   double Kb = 2*vMobT(B)*K(B)*vMobT(I)*K(I)/(vMobT(B)*K(B) + vMobT(I)*K(I));
00167   Point3D DX = m_mesh.getDX();
00168   double Pant = vPprev(I);
00169   double result;
00170   result=  por(I)*Bg(I)*P(I)/dt +
00171     (Kr*(P(I)-P(R)) + Kl*(P(I)-P(L)))/(DX(0)*DX(0)) +
00172     (Ku*(P(I)-P(U)) + Kd*(P(I)-P(D)))/(DX(1)*DX(1)) +
00173     (Kf*(P(I)-P(F)) + Kb*(P(I)-P(B)))/(DX(2)*DX(2))
00174     -
00175     (por(I)*Bg(I)*Pant/dt + por(I)*(St(I)-1.0)/dt);
00176     
00177   printf("Residual = %g\n", result);
00178     
00179 
00180   
00181 }