MixedHybridCompressibleNewton Class Reference

#include <mixedhybridcompressiblenewton.h>

Inheritance diagram for MixedHybridCompressibleNewton:

Collaboration diagram for MixedHybridCompressibleNewton:

Public Member Functions
	MixedHybridCompressibleNewton (OrthoMesh &mesh, Function3D &fPrescribedVelocities, Function3D &fPrescribedPression, const VecDouble &K, Function1D &density, LinearSolver &solver, int debugLevel)
	~MixedHybridCompressibleNewton ()
virtual void	iterate (TransportBase &trans)
Private Attributes
OrthoMesh &	m_mesh
Function1D &	m_density
const VecDouble &	m_vK
double	m_dt
SparsityPattern	m_sparsePattern
SparseMatrix< double >	m_A
VecDouble	m_B
VecDouble	m_Pn
VecDouble	m_Pnn
LinearSolver &	m_solver
int	_debugLevel
HybridFaceBC	m_bc

Detailed Description

MixedHybridCompressibleNewton implements a monophasic compressible flow in porous media through Newton iterations assuming isothermal flow. The incompressible flow equatios are given by B(p)*dp/dt + div(g(p) Vdt) = q Vdt = -K Grad(P) Where B(p) = dg/dp g(p) = the density of the fluid K = intrinsic permeability

Definition at line 17 of file mixedhybridcompressiblenewton.h.

Constructor & Destructor Documentation

MixedHybridCompressibleNewton::MixedHybridCompressibleNewton	(	OrthoMesh &	mesh,
		Function3D &	fPrescribedVelocities,
		Function3D &	fPrescribedPression,
		const VecDouble &	K,
		Function1D &	density,
		LinearSolver &	solver,
		int	debugLevel
	)

MixedHybridCompressibleNewton::~MixedHybridCompressibleNewton ( )

Member Function Documentation

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

Implements DynamicBase.

Definition at line 22 of file mixedhybridcompressiblenewton.cpp.

00024 00025   A = 0; 00026   B = 0; 00027 00028 00029 00030 00031 00032 00033   { 00034 00035 00036 00037 00038 00039     { 00040 00041 00042 00043 00044 00045 00046 00047 00048 00049 00050 00051 00052 00053 00054 00055 00056 00057 00058 00059 00060 00061 00062     } 00063   } 00064 00065 00066 00067 00068     { 00069 00070 00071 00072 00073 00074 00075 00076       { 00077 00078 00079 00080 00081 00082       } 00083 00084       { 00085 00086 00087 00088 00089       } 00090     } 00091 00092 }

class="fragment">00023 { //Zero entries //Get the face iterator and the number of faces. OrthoMesh::Face_It face = mesh.begin_face(); unsigned nFaces = mesh.numFaces(); //For each face for (unsigned faceIndex=0;faceIndex < nFaces; faceIndex++,face++) unsigned c1,c2; face->getAdjCellIndices(c1,c2); //Make sure the face is not at boundary if (!face->at_boundary()) const double P1 = Pn(c1); const double P2 = Pn(c2); const double K1 = K(c1); const double K2 = K(c2); const double Kh1 = K1*g(P1); const double Kh2 = K2*g(P2); const double Mh = 2*Kh1*Kh2/(Kh1+Kh2); //Harmonic average const double Keff =Mh*Inv_DX2; const double vv = Keff*(P1-P2); const double den = (Kh1+Kh2)*(Kh1+Kh2); const double d1 = Kh2*Kh2*K1*dg(P1)/den*(P1-P2)*Inv_DX2; const double d2 = Kh1*Kh1*K2*dg(P2)/den*(P1-P2)*Inv_DX2; A.add(c1,c1, d1+Keff); A.add(c1,c2, d2-Keff); A.add(c2,c2, -d2+Keff); A.add(c2,c1, -d1-Keff); B(c1)-=vv; B(c2)+=vv; //Now process the boundary conditions HybridFaceBC::iterator bcEnd = bc.end(); for(HybridFaceBC::iterator bcIt = bc.begin();bcIt!=bcEnd;bcIt++) OrthoMesh::Face_It face = bcIt->face; assert(face->at_boundary()); unsigned c1,c2; face->getAdjCellIndices(c1,c2); if (bcIt->bcType == Dirichlet) double vv double c = (c1 != OrthoMesh::invalidIndex()) ? c1 : c2; double Keff = 2*K(c)*vMobT(c)*face->area()*face->areaPerCellVol(); A.add(c,c,Keff); B(c) += Keff*bcIt->value; else //Neumann condition if (c1 != OrthoMesh::invalidIndex()) B(c1) += -bcIt->value; else B(c2) += +bcIt->value;