MixedHybridCompositionalSimplePw Class Reference

#include <mixedhybridcompositionalsimplepw.h>

Inheritance diagram for MixedHybridCompositionalSimplePw:

Collaboration diagram for MixedHybridCompositionalSimplePw:

Public Member Functions
	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)
virtual	~MixedHybridCompositionalSimplePw ()
virtual void	iterate (TransportBase &trans)
void	setScalarBC (OrthoMesh &mesh, ScalarBC &bc, Function3D &fPressure, Function3D &fTransBC, GeneralFunctionInterface &fPc, FlashCompositional &flash)
virtual void	rollback ()
Protected Member Functions
void	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)
Protected Attributes
VecDouble	m_vMobT
VecDouble	m_vGravTerm
VecDouble	m_vPc
VecDouble	m_vKPc
bool	m_bFirstIteration
ScalarBC	m_pcbc
GeneralFunctionInterface &	m_fpc

Detailed Description

Definition at line 11 of file mixedhybridcompositionalsimplepw.h.

Constructor & Destructor Documentation

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
	)

Definition at line 6 of file mixedhybridcompositionalsimplepw.cpp.

00007  :MixedHybridCompositional(mesh,fInitP,fPrescribedVelocities, fPrescribedPressure, K, vPor,fMobT,fFracFlux,grav,  dt,flash,solver,debugLevel, false),m_fpc(pc)
00008 {
00009   m_vGravTerm.reinit(m_mesh.numCells());
00010   m_vMobT.reinit(m_mesh.numCells());
00011   m_vPc.reinit(m_mesh.numCells());
00012   m_vKPc.reinit(m_mesh.numCells());
00013   m_bFirstIteration=true;
00014   m_solAnt=m_sol;
00015   setScalarBC(mesh,m_pcbc,fPrescribedPressure,fTransportBC,pc,flash);
00016   
00017 }

virtual MixedHybridCompositionalSimplePw::~MixedHybridCompositionalSimplePw ( ) [inline, virtual]

Definition at line 29 of file mixedhybridcompositionalsimplepw.h.

00029 {}

Member Function Documentation

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

Definition at line 23 of file mixedhybridcompositionalsimplepw.cpp.

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

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

Implements DynamicBase.

Definition at line 205 of file mixedhybridcompositionalsimplepw.cpp.

00206 {
00207   m_solAnt=m_sol;
00208   this->assemblySystem(m_A,m_B,m_dt,m_sol,m_vK,m_vPor,m_flash,m_fMobT,m_fFracFlux,m_fpc,m_grav);
00209   m_solver.solve(m_A,m_sol,m_B);
00210   //Point3D p(5,50,50);
00211   //OrthoMesh::Cell_It cell = m_mesh.getCellAt(p);
00212   //testSolution(cell,m_dt, v1,m_vPor,m_vK,m_sol,PAnt,m_vMobT,v2); 
00213   //Get the fluxes
00214   this->getNormalVelocitiesFromPressure(m_Vq,m_mesh,m_bc,m_vK,m_sol,m_vMobT,m_vGravTerm,m_pcbc,m_vKPc,m_vPc);
00215 
00216 }

void MixedHybridCompositionalSimplePw::rollback ( ) [virtual]

Reimplemented from CompressibleDynamic.

Definition at line 220 of file mixedhybridcompositionalsimplepw.cpp.

00221 {
00222   m_sol=m_solAnt;
00223 }

void MixedHybridCompositionalSimplePw::setScalarBC	(	OrthoMesh &	mesh,
		ScalarBC &	bc,
		Function3D &	fPressure,
		Function3D &	fTransBC,
		GeneralFunctionInterface &	fPc,
		FlashCompositional &	flash
	)

Definition at line 163 of file mixedhybridcompositionalsimplepw.cpp.

00164 {
00165   assert(flash.numPhases() == 2);
00166   assert(fTransBC.getImageDim() == flash.numComponents());
00167   unsigned ncmps = flash.numComponents();
00168   OrthoMesh::Face_It face = mesh.begin_face();
00169   OrthoMesh::Face_It endf = mesh.end_face();
00170   
00171   VecDouble cmps(flash.numComponents());
00172   VecDouble phasesVol(flash.numPhases());
00173   _ScalarBC pcbc;
00174   VecDouble vS(1);
00175   FlashData flashData(flash.numPhases(),flash.numComponents());
00176   for (;face!=endf;face++)
00177   {
00178     if (face->at_boundary())
00179     {
00180       Point3D p = face->barycenter();
00181       if (fTransBC.isInDomain(p))
00182       {
00183         if (fPressure.isInDomain(p))
00184         {
00185           for (unsigned i=0;i<ncmps;i++)
00186           {
00187             cmps(i)=fTransBC(p,i);
00188           }
00189           double P = fPressure(p);
00190           flash.flash(P,cmps,flashData);
00191           flash.getPhasesVolume(P,flashData,phasesVol);
00192           double factor=1.0/NumericMethods::vectorSum(phasesVol);
00193           cmps*=factor;
00194           vS(0)=phasesVol(0)/(phasesVol(0) + phasesVol(1));
00195           pcbc.value=fPc(vS);
00196           pcbc.face=face;
00197           bc.push_back(pcbc);
00198         }
00199       }
00200     }
00201   }
00202 }