MixedHybridIncompressible Class Reference

#include <mixedhybridincompressible.h>

Inheritance diagram for MixedHybridIncompressible:

Collaboration diagram for MixedHybridIncompressible:

Public Member Functions
	MixedHybridIncompressible (OrthoMesh &mesh, Function3D &fPrescribedVelocities, Function3D &fPrescribedPression, const VecDouble &K, GeneralFunctionInterface &fMobT, GeneralFunctionInterface &fFracFlux, double grav, VecDouble &densities, LinearSolver &solver, int debugLevel)
	~MixedHybridIncompressible ()
virtual void	iterate (TransportBase &trans)
virtual const VecDouble &	getPressureAtCells ()
virtual void	getNormalVelocityAtFaces (VecDouble &vFNC)
virtual void	printOutput ()
Protected Attributes
OrthoMesh &	m_mesh
GeneralFunctionInterface &	m_fMobT
GeneralFunctionInterface &	m_fFracFlux
VecDouble	m_densities
const VecDouble &	m_vK
VecDouble	m_vMobT
VecDouble	m_vGravTerm
double	m_grav
double	m_dt
SparsityPattern	m_sparsePattern
SparseMatrix< double >	m_A
VecDouble	m_B
VecDouble	m_sol
LinearSolver &	m_solver
int	_debugLevel
HybridFaceBC	m_bc

Detailed Description

Solve the pressure equation for incompressible flow. The equations are given by

div(Vdt) = q Vdt = - mobT(Sw) K (Grad(P) - (phow fracW(Sw) + phoO fracO(Sw)) grav _)

q = Source term mobT = Total Mobility phoW,phoO = fluid densities fracW,fracO = The fractional flux flows grav = Gravity vector

Definition at line 19 of file mixedhybridincompressible.h.

Constructor & Destructor Documentation

MixedHybridIncompressible::MixedHybridIncompressible	(	OrthoMesh &	mesh,
		Function3D &	fPrescribedVelocities,
		Function3D &	fPrescribedPression,
		const VecDouble &	K,
		GeneralFunctionInterface &	fMobT,
		GeneralFunctionInterface &	fFracFlux,
		double	grav,
		VecDouble &	densities,
		LinearSolver &	solver,
		int	debugLevel
	)

Definition at line 8 of file mixedhybridincompressible.cpp.

00009   :m_mesh(mesh),
00010    m_fMobT(fMobT),
00011    m_fFracFlux(fFracFlux),
00012    m_densities(densities),
00013    m_vK(K),
00014    m_grav(grav),
00015    m_solver(solver),
00016    _debugLevel(debugLevel)
00017 {
00018   
00019   
00020   this->setBoundaryCondition(m_bc,m_mesh,fPrescribedVelocities,fPrescribedPression,true);
00021   this->buildSparsityPattern(mesh,m_sparsePattern);
00022   m_A.reinit(m_sparsePattern);
00023   m_B.reinit(m_mesh.numCells());
00024   m_sol.reinit(m_mesh.numCells(),NAN);
00025   //  m_Vq.reinit(m_mesh.numFaces());
00026   m_vGravTerm.reinit(m_mesh.numCells());
00027   m_vMobT.reinit(m_mesh.numCells());
00028   
00029 
00030 }

MixedHybridIncompressible::~MixedHybridIncompressible ( ) [inline]

Definition at line 58 of file mixedhybridincompressible.h.

00058 {}

Member Function Documentation

virtual void MixedHybridIncompressible::getNormalVelocityAtFaces ( VecDouble & vFNC ) [inline, virtual]

Reimplemented from DynamicBase.

Definition at line 61 of file mixedhybridincompressible.h.

00061                                                         {  this->getNormalVelocitiesFromPressure(vFNC,m_mesh,m_bc,m_vK,m_sol,m_vMobT,m_vGravTerm);
00062   }

virtual const VecDouble& MixedHybridIncompressible::getPressureAtCells ( ) [inline, virtual]

Reimplemented from DynamicBase.

Definition at line 60 of file mixedhybridincompressible.h.

00060 {return m_sol;}

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

Implements DynamicBase.

Definition at line 33 of file mixedhybridincompressible.cpp.

00034 {
00035   assert(m_fFracFlux.getImageDim() == 1);
00036   assert(m_fMobT.getDomainDim() == 1);
00037   assert(m_fMobT.getImageDim() == 1);
00038   assert(m_fFracFlux.getDomainDim() == 1);
00039   assert(m_grav >= 0);
00040   
00041   //Allocate the vector of saturations. Since the sum
00042   //must be one the saturation of the last phase is
00043   //a function of the others, so we just have numPhases()-1
00044   //Independent saturation variables
00045   static VecDouble vS(1);
00046 
00047 
00048 
00049 
00050   //Get the cell iterator and the number of cells.
00051   OrthoMesh::Cell_It cell = m_mesh.begin_cell();
00052   Index nCells = m_mesh.numCells();
00053 
00054   //Write Zeros in the system
00055   m_A=0.0;
00056   m_B=0.0;
00057 
00058   ArrayOfVecDouble &transSol = trans.getSolutionAtCells();
00059   //For each cell do
00060   for (Index cellIndex=0;cellIndex < nCells; cellIndex++,cell++)
00061   {
00062     vS(0)=transSol(cellIndex,0);
00063     //transSol.getVecValues(cellIndex,&vS);
00064 
00065     m_vMobT(cellIndex)=m_fMobT(vS);
00066     
00067     
00068 
00069     double dMob = m_fFracFlux(vS,0);
00070     double mobAcum=dMob;
00071     double dGravTerm;
00072     dGravTerm=dMob*m_densities(0);
00073     for (unsigned k=1;k<vS.size();k++)
00074     {
00075       dMob = m_fFracFlux(vS,k);
00076       dGravTerm+=dMob*m_densities(k);
00077       mobAcum+=dMob;  
00078     }
00079     dGravTerm+=(1.0-mobAcum)*m_densities(vS.size());
00080 
00081     //Remember that the gravity is in the opposite direction
00082     //of the Z axis (grav acceleration vector  = <0,0,-grav>)
00083     //That is why we use the minus sign  here
00084     m_vGravTerm(cellIndex)=-dGravTerm*m_grav;
00085 
00086     assert(mobAcum<=1.0000000001);
00087   }
00088   MixedHybridBase::assemblySystem(m_A,m_B,m_mesh,m_bc,m_vK,m_vMobT,m_vGravTerm);
00089 
00090   //Solve the system, storing the solution on m_sol.
00091   m_solver.solve(m_A,m_sol,m_B);
00092 
00093   
00094   //Update the fluxes
00095   //this->getNormalVelocitiesFromPressure(m_Vq,m_mesh,m_bc,m_vK,m_sol,m_vMobT,m_vGravTerm);
00096 
00097 }

void MixedHybridIncompressible::printOutput ( ) [virtual]

Implements DynamicBase.

Definition at line 100 of file mixedhybridincompressible.cpp.

00101 {
00102   if (_debugLevel > 0)
00103   {
00104     VecDouble v(m_mesh.numCells());
00105     VecDouble Vv(m_mesh.numVertices());
00106     HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 
00107 
00108     m_mesh.projectCentralValuesAtVertices(getPressureAtCells(),Vv);
00109     hdf5.writeScalarField(Vv,"P");
00110     if (_debugLevel >= 3)
00111     {
00112       VecDouble vF(m_mesh.numFaces());
00113       getNormalVelocityAtFaces(vF);
00114       Matrix M(m_mesh.numVertices(),3);
00115       m_mesh.projectVelocitiesInFacesToVertices(M,vF);
00116       M.copyColumn(Vv,2);
00117       hdf5.writeScalarField(Vv,"VelZ");
00118       M.copyColumn(Vv,0);
00119       hdf5.writeScalarField(Vv,"VelX");
00120       M.copyColumn(Vv,1);
00121       hdf5.writeScalarField(Vv,"VelY");
00122 
00123       
00124     }
00125 
00126   }
00127   
00128 }