DoublePorosityDiff Class Reference

#include <doubleporositydiff.h>

Inheritance diagram for DoublePorosityDiff:
Inheritance graph
[legend]
Collaboration diagram for DoublePorosityDiff:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 DoublePorosityDiff (OrthoMesh &mesh, const VecDouble &K, const VecDouble &por, Function1D &fmob, Function1D &dPC, Function1D &_PC, BlockMatrixModule &blockMatrix)
virtual ~DoublePorosityDiff ()
virtual void iterate (double dt)
void iterate_linear (double dt)
double fracturePc (double sat)
VecDouble & fracturePot (const VecDouble &sat, VecDouble &C)
virtual void printOutput ()

Private Attributes

OrthoMesh & mesh
const VecDouble & _por
const VecDouble & _K
Function1D & _fmob
Function1D & _dPC
Function1D & _PC
VecDouble S_new
VecDouble _S_oldTime
VecDouble _vPot
BlockMatrixModule & _blockMatrix
double q_mf
double dq_mf
unsigned _numTimeStep

Detailed Description

Definition at line 11 of file doubleporositydiff.h.

Constructor & Destructor Documentation

DoublePorosityDiff::DoublePorosityDiff ( OrthoMesh &  mesh,
const VecDouble &  K,
const VecDouble &  por,
Function1D &  fmob,
Function1D &  dPC,
Function1D &  _PC,
BlockMatrixModule &  blockMatrix 
)

Definition at line 8 of file doubleporositydiff.cpp.

00009   :mesh(mesh),_por(por),_K(K),_fmob(fmob),_dPC(dPC),_PC(PC),_blockMatrix(blockMatrix)
00010 {
00011   _S_oldTime.reinit(mesh.numCells(),0);
00012   _S_oldTime= 0.023883;  // this value is the initial sturation (for 20% of water in the reservoir)
00013   q_mf = 0.0;
00014   dq_mf = 0.0;
00015   _numTimeStep=0;
00016 }

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

Definition at line 31 of file doubleporositydiff.h.

00031 {printf("I HAVE BEEN CALLED \n");}

Member Function Documentation

double DoublePorosityDiff::fracturePc ( double  sat  ) 

Definition at line 219 of file doubleporositydiff.cpp.

00220 {
00221   double gamma=20000, theta=100;                    //TODO, defined these as global var
00222   
00223   if (sat<=0)
00224     return 1; //TODO
00225   else if (sat>1)
00226     return 0;
00227   else
00228     return (1-sat)*(gamma/sat-gamma+theta);
00229 }

VecDouble & DoublePorosityDiff::fracturePot ( const VecDouble &  sat,
VecDouble &  C 
)

Definition at line 232 of file doubleporositydiff.cpp.

00233 {
00234   double rho_o = 1.8, rho_w=1, gravity=9.8, depht=0; //todo defined these as global var
00235   _vPot=C;
00236   for (int u=0;u<sat.size();u++)
00237   {
00238     _vPot(u)=_PC(sat(u))-(rho_o-rho_w)*gravity*depht;
00239   }
00240   return _vPot; 
00241 }

void DoublePorosityDiff::iterate ( double  dt  )  [virtual]

Implements DiffusiveStep.

Definition at line 19 of file doubleporositydiff.cpp.

00020 {
00021   // // TEST OF THE LINEAR ITERATION
00022   // // Comment the following function to run the newton's iteration
00023   // iterate_linear(dt); //solving the diffusive step without the nonlinearity-without the newton's iteration
00024   // printf("WARNING: THIS IS NOT THE NEWTON'S ITERATION \n");
00025   // return;
00026 
00027   //_S_oldTime.reinit(mesh.numCells(),0);
00028   // printf("S_oldTime=%g S_new=%g S_size=%u numCells=%u\n",_S_oldTime(9),S_new(9),_S_oldTime.size(),mesh.numCells());
00029   // return;
00030 
00031   //monitoring the number of time steps for diffusive step
00032   _numTimeStep=_numTimeStep+1;
00033   printf("Number of time steps for diffusive step= %u\n", _numTimeStep);
00034   //printf("monPC=%g\n",_PC(0.023883));
00035 
00036   //BEGINING OF THE NEWTON'S ITERATION
00037   VecDouble deltaS_oldIt(mesh.numCells());
00038   deltaS_oldIt=0.0;
00039   VecDouble deltaS_new(mesh.numCells());
00040   
00041   S_new.setRef(&(this->getTransport().getSolutionAtCells()(0,0)),mesh.numCells());
00042   
00043   VecDouble S_oldIt(mesh.numCells());               //saturation at the old iteration level
00044   S_oldIt = S_new;
00045 
00046   VecDouble S_transp(mesh.numCells());             //saturation from the transport level
00047   S_transp=S_new;
00048 
00049   //double q_mf=0.0, dq_mf=0.0;                                 //matrix-block term, not available
00050   double dS=0.001; //for the change in saturation
00051   
00052   double S_tol=1.e-3, dS_tol=1.e-3, S_error=0.0, deltaS_error=0.0;  // tolerance and error of convergence
00053   unsigned numIt=0;   //number of iteration
00054 
00055   VecDouble tmp(mesh.numCells());
00056   
00057   //constuction of the basis function from the matrix block
00058   for (unsigned u=0;u<_S_oldTime.size();u++)
00059   {
00060     // printf("S_old=%g PC(sold)=%g \n",_S_oldTime(u), _PC(_S_oldTime(u)));
00061     _S_oldTime(u) =_PC(_S_oldTime(u));//-(rho_o-rho_w)*gravity*depht;
00062   }
00063   _blockMatrix.construct_bases(_S_oldTime,dt);
00064   
00065   do 
00066   {
00067     S_error=0.0;
00068     deltaS_error=0.0;
00069   
00070     // loop to cells  
00071     for (OrthoMesh::Cell_It cell=mesh.begin_cell(); cell!=mesh.end_cell(); cell++)
00072     {
00073       double H=pow(cell->volume(),1/3.0);
00074       double sumXIsb=0.0,sumVb=0.0;
00075       for (unsigned neigID=0;neigID<6;neigID++) // loop over faces (or the neighboring cells)
00076       {
00077         if(cell->neighbor_index(neigID)==-1) //for boundary purposes, if we are at the boundary the flux is zero
00078         {
00079           continue;      
00080         }
00081         double Lsb=0.0,DiffCoef,effDiffCoef_F,KIsb, XIsb, effK; //saturation lagrange multipier at faces, intermediate coefs
00082                                                        //effective difusion ceofficient (at cell faces), effective permeability
00083         //saturation at faces
00084         Lsb=(_K(cell->index())*S_transp(cell->index())+_K(cell->neighbor_index(neigID))*S_transp(cell->neighbor_index(neigID)))/(_K(cell->index())+_K(cell->neighbor_index(neigID)));
00085         //effective permeability  
00086         effK=2*_K(cell->index())*_K(cell->neighbor_index(neigID))/(_K(cell->index())+_K(cell->neighbor_index(neigID)));
00087         effDiffCoef_F=effK*_fmob(Lsb)*_dPC(Lsb);
00088         DiffCoef=_K(cell->index())*_fmob(Lsb)*_dPC(Lsb); //diffusive coef
00089         //printf("effDiff=%g \n",effDiffCoef_F);
00090         assert(effDiffCoef_F <=0); //make sure that the diffusion coefficient is a negative value
00091         KIsb=1.0e-15; //-H/effDiffCoef_F; //TODO 
00092         //intermidiate coefficient        
00093         XIsb=2*DiffCoef/H;       
00094           
00095         sumXIsb = sumXIsb + XIsb/(1-KIsb*XIsb);
00096         sumVb = sumVb +  XIsb/(1-KIsb*XIsb)*(S_new(cell->index())-S_oldIt(cell->index())-deltaS_oldIt(cell->index())+(_K(cell->neighbor_index(neigID))/(_K(cell->index())+_K(cell->neighbor_index(neigID)))-KIsb*effDiffCoef_F/H)*(S_oldIt(cell->index())+deltaS_oldIt(cell->index())-S_oldIt(cell->neighbor_index(neigID))-deltaS_oldIt(cell->neighbor_index(neigID)) ));
00097       } //end of face loop
00098       // printf("S_oldIt %g S_new=%g \n",S_oldIt(cell->index()),S_new(cell->index()));
00099       //printf("PC(S_old)=%g sourceterm=%g der_source=%g \n",_S_oldTime(cell->index()),q_mf,dq_mf);
00100       deltaS_new(cell->index())=pow(_por(cell->index())/dt-sumXIsb/H-dq_mf,-1.0)*(-_por(cell->index())*(S_new(cell->index())-S_transp(cell->index()))/dt+sumVb/H + q_mf);
00101       //sum the saturation and delta_S for the update of the saturation
00102       S_new(cell->index())=S_new(cell->index())+deltaS_new(cell->index());
00103       S_error = S_error + pow(S_new(cell->index())-S_oldIt(cell->index()),2.0);
00104       deltaS_error = deltaS_error + pow(deltaS_new(cell->index()),2.0);
00105       //get the source term from the matrix block
00106       q_mf=_blockMatrix.get_source(cell->index(),_PC(S_new(cell->index())));
00107       //dq_mf =(_blockMatrix.get_source(cell->index(),_PC(S_new(cell->index()+dS)))-q_mf)/dS;
00108       // printf("PC(S_old)=%g PC(S_new)=%g sourceterm=%.10g der_source=%g \n",_S_oldTime(cell->index()),_PC(S_new(cell->index())),q_mf,dq_mf);
00109       // printf("S_oldIt %g S_new=%g \n",S_oldIt(cell->index()),S_new(cell->index()));
00110       //printf("S_new=%g \n",S_new(cell->index()));
00111       //assert(S_new(cell->index())<=1.1 && S_new(cell->index())>=-0.01);
00112     }//end of the cell loop
00113     
00114     S_error = pow(S_error,0.5);     //error estimation
00115     deltaS_error =  pow(deltaS_error,0.5);
00116     deltaS_oldIt = deltaS_new;
00117     S_oldIt = S_new;                 //set new values to old for new iteration if no convergence
00118 
00119     numIt = numIt + 1;              //count number of iteration 
00120     printf("Iteration=%u S_error=%.10g deltaS_error=%.10g \n",numIt, S_error,deltaS_error); //print number of iteration and the corresponding error
00121     
00122   } while(S_error > S_tol || deltaS_error >dS_tol);//end of convergence loop
00123   printf("END OF THE NEWTON'S ITERATE: THE SATURATION CONVERGES \n");
00124   // S_new.print(std::cout);
00125   _S_oldTime = S_new;
00126 }//end of funtion loop

void DoublePorosityDiff::iterate_linear ( double  dt  ) 

Definition at line 130 of file doubleporositydiff.cpp.

00131 {
00132   // FPcSquare mpc(300,0.2,0.15,1.0e+15);
00133   // printf("mpc026=%g\n",mpc(0.261306));
00134   // return;
00135   
00136   //monitoring the number of time steps for diffusive step
00137   _numTimeStep=_numTimeStep+1;
00138   printf("Number of time steps for diffusive step= %u\n", _numTimeStep);
00139   
00140   S_new.setRef(&(this->getTransport().getSolutionAtCells()(0,0)),mesh.numCells());
00141   
00142   VecDouble S_oldIt(mesh.numCells());               //saturation at the old iteration level
00143   S_oldIt = S_new;
00144 
00145 
00146   VecDouble S_transp(mesh.numCells());             //saturation from the transport level
00147   S_transp=S_new;
00148 
00149   //double q_mf=0.0;                                 //matrix-block term, not available
00150   
00151   double S_tol=1.e-3, S_error=0;  // tolerance and error of convergence
00152   unsigned numIt=0;   //number of iteration
00153 
00154   do 
00155   {
00156     S_error=0;
00157     //go over all the cells  
00158     for (OrthoMesh::Cell_It cell=mesh.begin_cell(); cell!=mesh.end_cell(); cell++)
00159     {
00160       double H=pow(cell->volume(),1/3.0);
00161       // printf("H=%g Cell_volume=%g \n",H,cell->volume());
00162       double sumXIsb=0.0,sumVb=0.0;
00163       for (unsigned neigID=0;neigID<6;neigID++) // loop over faces (or neighboring cells)
00164       {
00165         if(cell->neighbor_index(neigID)==-1) //for boundary purposes, if we are at the boundary the flux is zero
00166         {
00167           continue;      
00168         }
00169         double Lsb=0.0,DiffCoef,effDiffCoef_F,KIsb, XIsb, effK; //saturation lagrange multipier at faces, intermediate coefs
00170                                                        //effective difusion ceofficient (at cell faces), effective permeability
00171         //saturation at faces
00172         Lsb=(_K(cell->index())*S_transp(cell->index())+_K(cell->neighbor_index(neigID))*S_transp(cell->neighbor_index(neigID)))/(_K(cell->index())+_K(cell->neighbor_index(neigID)));
00173         //printf("S_transp=%g S_oldit=%g \n",S_transp(cell->index()),S_oldIt(cell->neighbor_index(neigID)));
00174         //effective permeability  
00175         effK=2*_K(cell->index())*_K(cell->neighbor_index(neigID))/(_K(cell->index())+_K(cell->neighbor_index(neigID)));
00176         //printf("effK=%g\n",effK);
00177         //effective diffusion coefficient
00178         //      printf("fmobility=%g devPC=%g \n",_fmob(Lsb),_dPC(Lsb));
00179         effDiffCoef_F=effK*_fmob(Lsb)*_dPC(Lsb);
00180         DiffCoef=_K(cell->index())*_fmob(Lsb)*_dPC(Lsb);
00181         //printf("effDiff=%g\n",effDiffCoef_F);
00182         //make sure that the diffusion coefficient is a negative value
00183         assert(effDiffCoef_F <=0);
00184         //diffusive coef        
00185         //
00186         KIsb=1.0e-5; //-H/effDiffCoef_F; //TODO
00187         //printf("KIsb=%.15g\n",KIsb);
00188         //intermidiate coefficient        
00189         XIsb=2*DiffCoef/H;
00190           
00191         sumXIsb = sumXIsb + XIsb/(1-KIsb*XIsb);
00192         sumVb = sumVb +  XIsb/(1-KIsb*XIsb)*(-S_oldIt(cell->index())+(_K(cell->neighbor_index(neigID))/(_K(cell->index())+_K(cell->neighbor_index(neigID)))-KIsb*effDiffCoef_F/H)*(S_oldIt(cell->index())-S_oldIt(cell->neighbor_index(neigID))));
00193       }//end face loop
00194             // //set the matrix source term
00195       // block_mat b_mat;
00196       // b_mat.basis_fcts(_S_oldTime(cell->index()));
00197       // //q_mf=b_mat.getSourceTerm(fracturePot(S_new(cell->index())));
00198       q_mf=0.0000;
00199       S_new(cell->index())=pow(_por(cell->index())/dt-sumXIsb/H,-1.0)*(_por(cell->index())*S_transp(cell->index())/dt+sumVb/H + q_mf);
00200       //printf("S_new=%g sourceterm=%g \n",S_new(cell->index()), q_mf);
00201       //assert(S_new(cell->index())<=1.1 && S_new(cell->index())>=-0.1);
00202       S_error = S_error + pow(S_new(cell->index())-S_oldIt(cell->index()),2.0);
00203       // printf("S_new=%g sourceterm=%g \n",S_new(cell->index()), q_mf);
00204     }//end of the cell loop
00205  
00206     S_error = pow(S_error,0.5);     //error estimation
00207     S_oldIt = S_new;                 //set new values to old for new iteration if no convergence
00208     numIt = numIt + 1;              //count number of iteration
00209     _S_oldTime = S_new;
00210     printf("Iteration=%u Error=%.10g \n",numIt, S_error); //print number of iteration and the corresponding error
00211     
00212   } while(S_error > S_tol);//end of convergence loop
00213   printf("WARNING: THIS IS NOT THE NEWTON'S ITERATION \n");
00214   printf("done, the saturation converges \n");
00215   // S_new.print(std::cout);
00216 }//end of funtion loop

void DoublePorosityDiff::printOutput (  )  [virtual]

Implements DiffusiveStep.

Definition at line 244 of file doubleporositydiff.cpp.

00245 {
00246   VecDouble vValues(mesh.numVertices());
00247   this->mesh.projectCentralValuesAtVertices(S_new,vValues);
00248   HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 
00249   hdf5.writeScalarField(vValues,"SDiff");
00250 }

Member Data Documentation

BlockMatrixModule& DoublePorosityDiff::_blockMatrix [private]

Definition at line 23 of file doubleporositydiff.h.

Function1D& DoublePorosityDiff::_dPC [private]

Definition at line 18 of file doubleporositydiff.h.

Function1D& DoublePorosityDiff::_fmob [private]

Definition at line 17 of file doubleporositydiff.h.

const VecDouble& DoublePorosityDiff::_K [private]

Definition at line 16 of file doubleporositydiff.h.

unsigned DoublePorosityDiff::_numTimeStep [private]

Definition at line 25 of file doubleporositydiff.h.

Function1D& DoublePorosityDiff::_PC [private]

Definition at line 19 of file doubleporositydiff.h.

const VecDouble& DoublePorosityDiff::_por [private]

Definition at line 15 of file doubleporositydiff.h.

VecDouble DoublePorosityDiff::_S_oldTime [private]

Definition at line 21 of file doubleporositydiff.h.

VecDouble DoublePorosityDiff::_vPot [private]

Definition at line 22 of file doubleporositydiff.h.

double DoublePorosityDiff::dq_mf [private]

Definition at line 24 of file doubleporositydiff.h.

OrthoMesh& DoublePorosityDiff::mesh [private]

Definition at line 14 of file doubleporositydiff.h.

double DoublePorosityDiff::q_mf [private]

Definition at line 24 of file doubleporositydiff.h.

VecDouble DoublePorosityDiff::S_new [private]

Definition at line 20 of file doubleporositydiff.h.

The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Sun Apr 8 23:13:02 2012 for CO2INJECTION by  doxygen 1.6.3