doubleporositydiff.cpp

Go to the documentation of this file.

#include "doubleporositydiff.h"
#include "hdf5orthowriter.h"

/*Domain Decomposition model (DDM) for diffusive step. The module solve the equation 55 in the paper "A parallelizable method for two-phase flows in naturally-fractured reservoirs: Computational Geosciences 1 (1997) 333–368" */


  
DoublePorosityDiff::DoublePorosityDiff(OrthoMesh &mesh,const VecDouble &K,const VecDouble &por, Function1D &fmob, Function1D & dPC,Function1D &PC, BlockMatrixModule& blockMatrix) 
  :mesh(mesh),_por(por),_K(K),_fmob(fmob),_dPC(dPC),_PC(PC),_blockMatrix(blockMatrix)
{
  _S_oldTime.reinit(mesh.numCells(),0);
  _S_oldTime= 0.023883;  // this value is the initial sturation (for 20% of water in the reservoir)
  q_mf = 0.0;
  dq_mf = 0.0;
  _numTimeStep=0;
}


void DoublePorosityDiff::iterate(double dt) 
{
  // // TEST OF THE LINEAR ITERATION
  // // Comment the following function to run the newton's iteration
  // iterate_linear(dt); //solving the diffusive step without the nonlinearity-without the newton's iteration
  // printf("WARNING: THIS IS NOT THE NEWTON'S ITERATION \n");
  // return;

  //_S_oldTime.reinit(mesh.numCells(),0);
  // printf("S_oldTime=%g S_new=%g S_size=%u numCells=%u\n",_S_oldTime(9),S_new(9),_S_oldTime.size(),mesh.numCells());
  // return;

  //monitoring the number of time steps for diffusive step
  _numTimeStep=_numTimeStep+1;
  printf("Number of time steps for diffusive step= %u\n", _numTimeStep);
  //printf("monPC=%g\n",_PC(0.023883));

  //BEGINING OF THE NEWTON'S ITERATION
  VecDouble deltaS_oldIt(mesh.numCells());
  deltaS_oldIt=0.0;
  VecDouble deltaS_new(mesh.numCells());
  
  S_new.setRef(&(this->getTransport().getSolutionAtCells()(0,0)),mesh.numCells());
  
  VecDouble S_oldIt(mesh.numCells());               //saturation at the old iteration level
  S_oldIt = S_new;

  VecDouble S_transp(mesh.numCells());             //saturation from the transport level
  S_transp=S_new;

  //double q_mf=0.0, dq_mf=0.0;                                 //matrix-block term, not available
  double dS=0.001; //for the change in saturation
  
  double S_tol=1.e-3, dS_tol=1.e-3, S_error=0.0, deltaS_error=0.0;  // tolerance and error of convergence
  unsigned numIt=0;   //number of iteration

  VecDouble tmp(mesh.numCells());
  
  //constuction of the basis function from the matrix block
  for (unsigned u=0;u<_S_oldTime.size();u++)
  {
    // printf("S_old=%g PC(sold)=%g \n",_S_oldTime(u), _PC(_S_oldTime(u)));
    _S_oldTime(u) =_PC(_S_oldTime(u));//-(rho_o-rho_w)*gravity*depht;
  }
  _blockMatrix.construct_bases(_S_oldTime,dt);
  
  do 
  {
    S_error=0.0;
    deltaS_error=0.0;
  
    // loop to cells  
    for (OrthoMesh::Cell_It cell=mesh.begin_cell(); cell!=mesh.end_cell(); cell++)
    {
      double H=pow(cell->volume(),1/3.0);
      double sumXIsb=0.0,sumVb=0.0;
      for (unsigned neigID=0;neigID<6;neigID++) // loop over faces (or the neighboring cells)
      {
        if(cell->neighbor_index(neigID)==-1) //for boundary purposes, if we are at the boundary the flux is zero
        {
          continue;      
        }
        double Lsb=0.0,DiffCoef,effDiffCoef_F,KIsb, XIsb, effK; //saturation lagrange multipier at faces, intermediate coefs
                                                       //effective difusion ceofficient (at cell faces), effective permeability
        //saturation at faces
        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)));
        //effective permeability  
        effK=2*_K(cell->index())*_K(cell->neighbor_index(neigID))/(_K(cell->index())+_K(cell->neighbor_index(neigID)));
        effDiffCoef_F=effK*_fmob(Lsb)*_dPC(Lsb);
        DiffCoef=_K(cell->index())*_fmob(Lsb)*_dPC(Lsb); //diffusive coef
        //printf("effDiff=%g \n",effDiffCoef_F);
        assert(effDiffCoef_F <=0); //make sure that the diffusion coefficient is a negative value
        KIsb=1.0e-15; //-H/effDiffCoef_F; //TODO 
        //intermidiate coefficient        
        XIsb=2*DiffCoef/H;       
          
        sumXIsb = sumXIsb + XIsb/(1-KIsb*XIsb);
        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)) ));
      } //end of face loop
      // printf("S_oldIt %g S_new=%g \n",S_oldIt(cell->index()),S_new(cell->index()));
      //printf("PC(S_old)=%g sourceterm=%g der_source=%g \n",_S_oldTime(cell->index()),q_mf,dq_mf);
      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);
      //sum the saturation and delta_S for the update of the saturation
      S_new(cell->index())=S_new(cell->index())+deltaS_new(cell->index());
      S_error = S_error + pow(S_new(cell->index())-S_oldIt(cell->index()),2.0);
      deltaS_error = deltaS_error + pow(deltaS_new(cell->index()),2.0);
      //get the source term from the matrix block
      q_mf=_blockMatrix.get_source(cell->index(),_PC(S_new(cell->index())));
      //dq_mf =(_blockMatrix.get_source(cell->index(),_PC(S_new(cell->index()+dS)))-q_mf)/dS;
      // 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);
      // printf("S_oldIt %g S_new=%g \n",S_oldIt(cell->index()),S_new(cell->index()));
      //printf("S_new=%g \n",S_new(cell->index()));
      //assert(S_new(cell->index())<=1.1 && S_new(cell->index())>=-0.01);
    }//end of the cell loop
    
    S_error = pow(S_error,0.5);     //error estimation
    deltaS_error =  pow(deltaS_error,0.5);
    deltaS_oldIt = deltaS_new;
    S_oldIt = S_new;                 //set new values to old for new iteration if no convergence

    numIt = numIt + 1;              //count number of iteration 
    printf("Iteration=%u S_error=%.10g deltaS_error=%.10g \n",numIt, S_error,deltaS_error); //print number of iteration and the corresponding error
    
  } while(S_error > S_tol || deltaS_error >dS_tol);//end of convergence loop
  printf("END OF THE NEWTON'S ITERATE: THE SATURATION CONVERGES \n");
  // S_new.print(std::cout);
  _S_oldTime = S_new;
}//end of funtion loop

#include "fpcsquare.h"

void DoublePorosityDiff::iterate_linear(double dt) 
{
  // FPcSquare mpc(300,0.2,0.15,1.0e+15);
  // printf("mpc026=%g\n",mpc(0.261306));
  // return;
  
  //monitoring the number of time steps for diffusive step
  _numTimeStep=_numTimeStep+1;
  printf("Number of time steps for diffusive step= %u\n", _numTimeStep);
  
  S_new.setRef(&(this->getTransport().getSolutionAtCells()(0,0)),mesh.numCells());
  
  VecDouble S_oldIt(mesh.numCells());               //saturation at the old iteration level
  S_oldIt = S_new;


  VecDouble S_transp(mesh.numCells());             //saturation from the transport level
  S_transp=S_new;

  //double q_mf=0.0;                                 //matrix-block term, not available
  
  double S_tol=1.e-3, S_error=0;  // tolerance and error of convergence
  unsigned numIt=0;   //number of iteration

  do 
  {
    S_error=0;
    //go over all the cells  
    for (OrthoMesh::Cell_It cell=mesh.begin_cell(); cell!=mesh.end_cell(); cell++)
    {
      double H=pow(cell->volume(),1/3.0);
      // printf("H=%g Cell_volume=%g \n",H,cell->volume());
      double sumXIsb=0.0,sumVb=0.0;
      for (unsigned neigID=0;neigID<6;neigID++) // loop over faces (or neighboring cells)
      {
        if(cell->neighbor_index(neigID)==-1) //for boundary purposes, if we are at the boundary the flux is zero
        {
          continue;      
        }
        double Lsb=0.0,DiffCoef,effDiffCoef_F,KIsb, XIsb, effK; //saturation lagrange multipier at faces, intermediate coefs
                                                       //effective difusion ceofficient (at cell faces), effective permeability
        //saturation at faces
        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)));
        //printf("S_transp=%g S_oldit=%g \n",S_transp(cell->index()),S_oldIt(cell->neighbor_index(neigID)));
        //effective permeability  
        effK=2*_K(cell->index())*_K(cell->neighbor_index(neigID))/(_K(cell->index())+_K(cell->neighbor_index(neigID)));
        //printf("effK=%g\n",effK);
        //effective diffusion coefficient
        //      printf("fmobility=%g devPC=%g \n",_fmob(Lsb),_dPC(Lsb));
        effDiffCoef_F=effK*_fmob(Lsb)*_dPC(Lsb);
        DiffCoef=_K(cell->index())*_fmob(Lsb)*_dPC(Lsb);
        //printf("effDiff=%g\n",effDiffCoef_F);
        //make sure that the diffusion coefficient is a negative value
        assert(effDiffCoef_F <=0);
        //diffusive coef        
        //
        KIsb=1.0e-5; //-H/effDiffCoef_F; //TODO
        //printf("KIsb=%.15g\n",KIsb);
        //intermidiate coefficient        
        XIsb=2*DiffCoef/H;
          
        sumXIsb = sumXIsb + XIsb/(1-KIsb*XIsb);
        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))));
      }//end face loop
            // //set the matrix source term
      // block_mat b_mat;
      // b_mat.basis_fcts(_S_oldTime(cell->index()));
      // //q_mf=b_mat.getSourceTerm(fracturePot(S_new(cell->index())));
      q_mf=0.0000;
      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);
      //printf("S_new=%g sourceterm=%g \n",S_new(cell->index()), q_mf);
      //assert(S_new(cell->index())<=1.1 && S_new(cell->index())>=-0.1);
      S_error = S_error + pow(S_new(cell->index())-S_oldIt(cell->index()),2.0);
      // printf("S_new=%g sourceterm=%g \n",S_new(cell->index()), q_mf);
    }//end of the cell loop
 
    S_error = pow(S_error,0.5);     //error estimation
    S_oldIt = S_new;                 //set new values to old for new iteration if no convergence
    numIt = numIt + 1;              //count number of iteration
    _S_oldTime = S_new;
    printf("Iteration=%u Error=%.10g \n",numIt, S_error); //print number of iteration and the corresponding error
    
  } while(S_error > S_tol);//end of convergence loop
  printf("WARNING: THIS IS NOT THE NEWTON'S ITERATION \n");
  printf("done, the saturation converges \n");
  // S_new.print(std::cout);
}//end of funtion loop


double DoublePorosityDiff::fracturePc(double sat)
{
  double gamma=20000, theta=100;                    //TODO, defined these as global var
  
  if (sat<=0)
    return 1; //TODO
  else if (sat>1)
    return 0;
  else
    return (1-sat)*(gamma/sat-gamma+theta);
}


  VecDouble& DoublePorosityDiff::fracturePot(const VecDouble &sat,VecDouble &C)
{
  double rho_o = 1.8, rho_w=1, gravity=9.8, depht=0; //todo defined these as global var
  _vPot=C;
  for (int u=0;u<sat.size();u++)
  {
    _vPot(u)=_PC(sat(u))-(rho_o-rho_w)*gravity*depht;
  }
  return _vPot; 
}


void DoublePorosityDiff::printOutput() 
{
  VecDouble vValues(mesh.numVertices());
  this->mesh.projectCentralValuesAtVertices(S_new,vValues);
  HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 
  hdf5.writeScalarField(vValues,"SDiff");
}

---