blockmatrixmodule.cpp

Go to the documentation of this file.

#include "blockmatrixmodule.h"
#include <stdio.h>
#include<fstream>
#include "umfpacksolver.h"
#include <iostream>
#include "hdf5orthowriter.h"
#include "xdmfwriter.h"

using namespace std;

using namespace dealii;

void BlockMatrixModule::iterate(TransportBase &trans) 
{
  //  n_fracture_cells=125 * 2;
  // cout<<"frac cell "<<n_fracture_cells<<endl;
  VecDouble P;
  P.reinit(n_fracture_cells);
  for (unsigned i=0; i < n_fracture_cells; i++)
    P(i)=0.2;

  
  construct_bases (P, _dt);
  
  // std::ofstream out ("sp1");
  // sp.print_gnuplot (out);
  // abort();
}

void BlockMatrixModule::buildSparsityPattern() 
{
  // SparsityPattern sp;
  // unsigned m=_mesh.numRawVertices();
  // sp.reinit(2*m,2*m,14);
 
  // printf(" Vertices: sparsity %u \n", 2*_m);

  OrthoMesh::Raw_Vertex_It end = _mesh.end_raw_vertice();

  unsigned ix, iy, iz, ver_xy, _ver_xy,ver_tol;
  _ver_xy=(_mesh.numElemX()+1)*(_mesh.numElemY()+1); // total # of vertices in xy-plane
  ver_xy= (_mesh.numElemX()+1)*(_mesh.numElemY()+1)-(_mesh.numElemX()+1); // number used for front and back face boundaries.
  // printf("xy %u \n", ver_xy);
  ver_tol = (_mesh.numElemX()+1)*(_mesh.numElemY()+1)*(_mesh.numElemY()+1);  

 
  /*Front face boundary*/
  OrthoMesh::Raw_Vertex_It vf = _mesh.begin_raw_vertice();

  for (iz=0; iz<_mesh.numElemZ()+1; iz++)
  {
    //  printf("check vertex %u \n ", vf->index());
    for (ix=0; ix<_mesh.numElemX()+1; ix++)
    {
      sp.add(vf->index(),vf->index());
      sp.add(_m+vf->index(),_m+vf->index());
      //      printf("check vertex %u \n ", vf->index());
      vf++;      
    }
    
    for (unsigned k=0; k < ver_xy; k++)
    {
          vf++;
          //  printf("check vertex %u \n ", vf->index());
    }
  }
 
  /*Back face boundary*/
  OrthoMesh::Raw_Vertex_It vb = _mesh.begin_raw_vertice();
  // printf("here %u \n", vb->index());
 
  for (iz=0; iz<_mesh.numElemZ()+1; iz++)
  {
     for (unsigned k=0; k < ver_xy; k++)
          vb++;
     //  printf("check vertex vb  %u \n ", vb->index());
    for (ix=0; ix<_mesh.numElemX()+1; ix++)
    {
      sp.add(vb->index(),vb->index());
      sp.add(_m+vb->index(),_m+vb->index());
      //  printf("check vertex vb  %u \n ", vb->index());
      vb++;     
    }
  }

  /*Left face boundary*/
  OrthoMesh::Raw_Vertex_It vl = _mesh.begin_raw_vertice();
   
  for (iz=0; iz<_mesh.numElemZ()+1; iz++)
  {    
    //  for (iy=0; iy<_mesh.numElemY()+1; iy++)
    for (iy=0; iy<_mesh.numElemY()+1; iy++)
    {  
      sp.add(vl->index(),vl->index());
      sp.add(_m+vl->index(),_m+vl->index());
      //  printf("vl %u \n", vl->index()); 
      for (unsigned k=0; k <_mesh.numElemX()+1; k++)
         vl++;
          
    }
  }

  /*Right face boundary*/  
  OrthoMesh::Raw_Vertex_It vr = _mesh.begin_raw_vertice();
  // printf("here vr %u \n", vr->index());
  for (unsigned k=0; k <_mesh.numElemX(); k++)
         vr++;
  //     printf("VR: out %u \n", vr->index());
     
  for (iz=0; iz<_mesh.numElemZ()+1; iz++)
  {
    for (iy=0; iy<_mesh.numElemY()+1; iy++)
    {    
      sp.add(vr->index(),vr->index());
      sp.add(_m+vr->index(),_m+vr->index());
      //       cout<<"vr \t"<< vr->index()<<endl;
      for (unsigned k=0; k <_mesh.numElemX()+1; k++)
      {
        
        if (vr->index() != _m-1)
        {
        vr++;
        //      cout<<"vertices inside\t"<< vr->index()<<endl;
        }
      }   
    }
  }


  /*Bottom face boundary*/ 
  OrthoMesh::Raw_Vertex_It vbt = _mesh.begin_raw_vertice();
 
  iz=0;
 
  for (iy=0; iy<_mesh.numElemY()+1; iy++)
  {    
    for (ix=0; ix<_mesh.numElemX()+1; ix++)
    {  
      sp.add(vbt->index(),vbt->index());
      sp.add(_m+vbt->index(),_m+vbt->index());
      //        printf("vbt %u \n", vbt->index());        
     vbt++;
    }
  }
 
  /*Top (up) face boundary*/  
  OrthoMesh::Raw_Vertex_It vu = _mesh.begin_raw_vertice();
  
  for (unsigned k=0 ; k <_m-(_ver_xy); k++)
    vu++;
  // printf("VU %u \n", vu->index());
  
  iz=_mesh.numElemZ()+1;
  //   if (I == -1)
  for (iy=0; iy<_mesh.numElemY()+1; iy++)
  {   
    for (ix=0; ix<_mesh.numElemX()+1; ix++)
    {
      sp.add(vu->index(),vu->index());
      sp.add(_m+vu->index(),_m+vu->index());
      //  printf("vu %u \n", vu->index());        
      vu++;
    }
  }

  /*Center vertices*/
  OrthoMesh::Raw_Vertex_It vc = _mesh.begin_raw_vertice();
  
  for (int i=0; i < (_mesh.numElemX()+1) * (_mesh.numElemY()+1) + (_mesh.numElemX()+2); i++)
  {
    vc++;
  }
  // printf("vc %u\n", vc->index());
  //--------------------------------------------------------------------
   for (iz=1; iz < _mesh.numElemZ(); iz++)
   {
    for (iy=1; iy < _mesh.numElemY(); iy++)
    {
      for (ix=1; ix < _mesh.numElemX(); ix++)
      {
        //      printf("vc %u\n", vc->index());
          //si_c 1st eq.
          sp.add(vc->index(),vc->index());
          //si_w 1st eq.
          sp.add(vc->index(),_m+vc->index());
          sp.add(vc->index(),_m+vc->index_left());
          sp.add(vc->index(), _m+ vc->index_right());
          sp.add(vc->index(),_m+vc->index_bottom());
          sp.add(vc->index(),_m+vc->index_up());
          sp.add(vc->index(),_m+ vc->index_back());
          sp.add(vc->index(),_m+ vc->index_front());
          //si_c 2nd eq.
          sp.add(_m+vc->index(),vc->index());
          sp.add(_m+vc->index(),vc->index_left());
          sp.add(_m+vc->index(), vc->index_right());
          sp.add(_m+vc->index(),vc->index_bottom());
          sp.add(_m+vc->index(),vc->index_up());
          sp.add(_m+vc->index(), vc->index_back());
          sp.add(_m+vc->index(), vc->index_front());
          //si_w 2nd eq.
          sp.add(_m+vc->index(),_m+vc->index());
          sp.add(_m+vc->index(),_m+vc->index_left());
          sp.add(_m+vc->index(),_m+ vc->index_right());
          sp.add(_m+vc->index(),_m+vc->index_bottom());
          sp.add(_m+vc->index(),_m+vc->index_up());
          sp.add(_m+vc->index(),_m+ vc->index_back());
          sp.add(_m+vc->index(),_m+ vc->index_front());       
          vc++;
        }

      for (unsigned i=0; i<2; i++)
        vc++;
    }

    for (unsigned i=0; i< 2*(_mesh.numElemX()+1); i++)
        vc++;
 
    } 
 



  //-----------------------------------------------------------------------

  /*  
   for (iz=0; iz < _mesh.numElemZ()+1; iz++)
    for (iy=0; iy < _mesh.numElemY()+1; iy++)
      for (ix=0; ix < _mesh.numElemX()+1; ix++)
      {
        //      cout<<"here\t"<<vc->index()<<endl;
        if (ix != 0 && iy != 0 && iz != 0 && ix != _mesh.numElemX() &&  iy != _mesh.numElemY() &&  iz != _mesh.numElemZ())
        {
          // printf("vc %u %u %u %u %u\n", vc->index(), vc->index_bottom(), ix, iy, iz);
          //si_c 1st eq.
          sp.add(vc->index(),vc->index());
          //si_w 1st eq.
          sp.add(vc->index(), _m+vc->index());
          sp.add(vc->index(), _m+vc->index_left());
          sp.add(vc->index(), _m+vc->index_right());
          sp.add(vc->index(), _m+vc->index_bottom());
          sp.add(vc->index(), _m+vc->index_up());
          sp.add(vc->index(), _m+vc->index_back());
          sp.add(vc->index(), _m+vc->index_front());
          //si_c 2nd eq.
          sp.add(_m+vc->index(),vc->index());
          sp.add(_m+vc->index(),vc->index_left());
          sp.add(_m+vc->index(), vc->index_right());
          sp.add(_m+vc->index(),vc->index_bottom());
          sp.add(_m+vc->index(),vc->index_up());
          sp.add(_m+vc->index(), vc->index_back());
          sp.add(_m+vc->index(), vc->index_front());
          //si_w 2nd eq.
          sp.add(_m+vc->index(),_m+vc->index());
          sp.add(_m+vc->index(),_m+vc->index_left());
          sp.add(_m+vc->index(),_m+ vc->index_right());
          sp.add(_m+vc->index(),_m+vc->index_bottom());
          sp.add(_m+vc->index(),_m+vc->index_up());
          sp.add(_m+vc->index(),_m+ vc->index_back());
          sp.add(_m+vc->index(),_m+ vc->index_front());
          printf("vc %u \n", vc->index());
        }
        vc++;
        
        
      }
  */
     sp.compress();
  
     //      printf("Sparsity end: %u \n", sp.n_nonzero_elements());
       
     // std::ofstream out ("sp1");
     //  sp.print_gnuplot (out);
     /* to save the image:
        set terminal png
        set output 'sp1.png'
        plot "sp1" with points */
           
}

void BlockMatrixModule::assembleSystem(const VecDouble &sat_old, const VecDouble &Sol_old,double _k, double _por, double _bc_pc, double dt) 
{
  
   OrthoMesh::Raw_Vertex_It vt=_mesh.begin_raw_vertice();
   OrthoMesh::Raw_Vertex_It vend=_mesh.end_raw_vertice();
   // SparseMatrix<double> M;
   // cout<<"Assemble System"<<endl;
      M.reinit(sp);

      //     rhs_Si.reinit(2*m);
    
           
      while (vt != vend)
      {
        
            OrthoVerticeAccessor ortho(_mesh,vt->index());
            if(ortho.at_boundary())
            {
              // printf("Boundary \n");
              M.add(vt->index(),vt->index(),1);
              M.add(_m+vt->index(),_m+vt->index(),1);
              // cout<<"index \t"<<vt->index()<<endl;
              // if (vt->index() <_m)
              //  {
                rhs_Si(vt->index())=_bc_pc; //Eq. 46 boundary condition
                rhs_Sihom(vt->index())=1;  //Eq. 47
                // }
              // else
              // {
                 rhs_Si(_m+vt->index())=0;
                 rhs_Sihom(_m+vt->index())=0;
                 //       }
                
            }
            else
            {
             //  printf("Not boundary \n");
             //  cout<<"here\t"<<vt->index()<<endl;
              //              Sol_old.print();
               M.add(vt->index(),vt->index(), _por /( _fdpc(sat_old(vt->index())) *_dt ) ); // si_c eq. 46
               //saturation at the midpoint
               // cout<<"sat old \t "<<vt->index()<<'\t'<<sat_old(vt->index())<<endl;
              
               double L = 0.5*( sat_old( vt->index()) + sat_old(vt->index_left()) );
               double R = 0.5*(  sat_old( vt->index()) + sat_old(vt->index_right()) );
               double F = 0.5*( sat_old( vt->index()) + sat_old(vt->index_front()) );
               double B = 0.5*( sat_old( vt->index()) + sat_old(vt->index_back()) );
               double U  = 0.5*(  sat_old( vt->index()) + sat_old(vt->index_up()) );
               double Bt = 0.5*( sat_old( vt->index()) + sat_old(vt->index_bottom()) );
               double _Dx = _dx * _dx;
               double _Dy = _dy * _dy;
               double _Dz = _dz * _dz;
               //              cout<<"L \t"<<L<<'\t'<<_fmobt(L)-_fmobrW(L)<<'\t'<<endl;
              
                M.add(vt->index(),_m+vt->index_left(), - _k *  _fmobrW(L)/_Dx );
                M.add(vt->index(),_m+vt->index_right(),-_k *  _fmobrW(R)/_Dx );
                M.add(vt->index(),_m+vt->index_front(),-_k *  _fmobrW(F)/_Dy );
                M.add(vt->index(),_m+vt->index_back(),- _k *  _fmobrW(B)/_Dy );
                M.add(vt->index(),_m+vt->index_bottom(),- _k * _fmobrW(Bt)/_Dz );
                M.add(vt->index(),_m+vt->index_up(), - _k * _fmobrW(U)/_Dz );
                M.add(vt->index(),_m+vt->index(), _k *(  (_fmobrW(L) + _fmobrW(R))/_Dx + (_fmobrW(F) + _fmobrW(B))/_Dy + (_fmobrW(Bt) + _fmobrW(U))/_Dz ) );
                rhs_Si(vt->index())=_por *Sol_old(vt->index()) /( _fdpc( sat_old(vt->index()) )*_dt);
                //cout<<"old solution" << sat_old(vt->index())<<endl;
                rhs_Sihom(vt->index())=0;
                //      ------------Eq. 47----------------------
                //      Si_c
                //      cout<<"ver \t"<<vt->index()<<'\t'<<_m+vt->index()<<'\t'<<_fmobrW(L) + _fmobrW(R) + _fmobrW(F) + _fmobrW(B) + _fmobrW(Bt) + _fmobrW(U)<<endl;
                M.add(_m+vt->index(),vt->index(), _k *( (_fmobt(L)-_fmobrW(L))/_Dx +  (_fmobt(R)-_fmobrW(R))/_Dx +  (_fmobt(F)-_fmobrW(F))/_Dy + (_fmobt(B)-_fmobrW(B))/_Dy +  (_fmobt(Bt)-_fmobrW(Bt))/_Dz + (_fmobt(U)-_fmobrW(U))/_Dz ) );
                M.add(_m+vt->index(),vt->index_left(),- _k * (_fmobt(L)- _fmobrW(L))/_Dx );
                M.add(_m+vt->index(),vt->index_right(),- _k  * (_fmobt(R)- _fmobrW(R))/_Dx );
                M.add(_m+vt->index(),vt->index_front(), -_k * (_fmobt(F)- _fmobrW(F))/_Dy );
                M.add(_m+vt->index(),vt->index_back(), - _k * (_fmobt(B)- _fmobrW(B))/_Dy );
                M.add(_m+vt->index(),vt->index_bottom(),- _k  * (_fmobt(Bt)- _fmobrW(Bt))/_Dz );
                M.add(_m+vt->index(),vt->index_up(), - _k * (_fmobt(U)- _fmobrW(U))/_Dz );

                //Si_w
                M.add(_m+vt->index(),_m+vt->index(), _k *( (_fmobt(L)+_fmobt(R))/_Dx + (_fmobt(F)+_fmobt(B))/_Dy + (_fmobt(Bt)+_fmobt(U))/_Dz ));
                M.add(_m+vt->index(),_m+vt->index_left(), -_k *_fmobt(L)/_Dx);
                M.add(_m+vt->index(),_m+vt->index_right(), -_k *_fmobt(R)/_Dx );
                M.add(_m+vt->index(),_m+vt->index_front(), -_k *_fmobt(F)/_Dy);
                M.add(_m+vt->index(),_m+vt->index_back(), -_k *_fmobt(B)/_Dy );
                M.add(_m+vt->index(),_m+vt->index_bottom(), -_k *_fmobt(Bt)/_Dz );
                M.add(_m+vt->index(),_m+vt->index_up(), -_k *_fmobt(U)/_Dz );
                
                
        
            } 
            vt++;
       }

      // FILE *fdata;
      // fdata= fopen("fdata","w");
      //    for (unsigned i=0; i< 2*_m; i++)
      //    {
      //           for (unsigned j=0; j< 2*_m; j++)
      //           {
      //             //printf("(%u,%u)%g \n",i,j, M.el(i,j));
             
      //              fprintf(fdata,"%g ",M.el(i,j));
      //              //    if (sp.exists(i,j))
      //              //        cout<<'\t'<<i<<','<<j<<'\t'<<M.operator()(i,j);
      //           }
      //           //      cout<<endl;
      //         }
      //         cout<<"total number of rows"<<_m<<endl;
      
      Sol_Si();
       Sol_Sihom();
       // rhs_Si.print();
       // rhs_Sihom.print();
           //     // cout<<endl;
       //   Solution1.print();
      //cout<<endl;
         //     Solution2.print();
      //    cout<<Solution1(40)<<endl;
       //   for (unsigned i=0; i< 2*_m ;i++)
  // {
  //    cout<<"rhs \t"<<i<<'\t'<<rhs_Sihom(i)<<endl;
  //       cout<<i<<'\t'<<Solution1(i)<<endl;
  // }            
}


void BlockMatrixModule::Sol_Si() 
{
  
  UMFPACKSolver umfp;
  umfp.solve(M,Solution1,rhs_Si);

  //  Solution
  for (unsigned i=0; i< _m ;i++)
  {
     Sol(i)=Solution1(i);
    //     cout<<Sol(i)<<endl;
  }

    
}

void BlockMatrixModule::Sol_Sihom() 
{
  UMFPACKSolver umfp;
  umfp.solve(M,Solution2,rhs_Sihom);
  for (unsigned i=0; i< _m ; i++)
    Sol_hom(i)=Solution2(i);

 }

void BlockMatrixModule::construct_bases(VecDouble &bc_pc,double dt) 
{
  //_dt=0.1;
  VecDouble SiPrev(_m);
  VecDouble SwPrev;
 
  //For each fracture block 
  for (unsigned blockId =0; blockId < n_fracture_cells ; blockId++)
  {
    //SiPrev = Si + (PcNew-PcOld)*SiHom
    SiPrev.copyRow(blockId,m_MSi); //prev capillary pressure solution
    //  cout<<"construct basis \n";
   
    VecDoubleRef SiHom(&(m_MHomSi(blockId,0)),_m);
  
    SiPrev.sadd(1,_BcPcNew(blockId)-_BcPcOld(blockId),SiHom);
  
    //Calculate the saturation corresponding to the SiPrev      
    SwPrev = SiPrev;
    _fpcinv.evaluate(SwPrev); // saturation at previous time level

    //  SwPrev.print();
    //Store the integral of the sw field
    _SwIntegral(blockId)=_mesh.getIntegralAtVertices(SwPrev);
    
    VecDoubleRef Si(&(m_MSi(blockId,0)),_m);

    create_bases_functions(bc_pc(blockId),SwPrev,SiPrev,_vPor(blockId),_vK(blockId),Si,SiHom);
    _BcPcOld(blockId)=bc_pc(blockId);

    HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter();
    hdf5.writeScalarField(Si,"Psi","MBmesh");

  }
    
    //  for (unsigned i=0; i<n_fracture_cells; i++)
    // {
    //   for (unsigned j=0; j<_m; j++)
    //     cout<<j<<'\t'<<m_MSi(i,j)<<endl;
    //   cout<<endl<<"new row"<<endl;
    // }

  // FILE *fout;
  //    fout= fopen("source","w");
           
  //      for (unsigned i=0; i< 200; i++)
  //      {
  //     fprintf(fout, "%u %g \n",i,get_source(0,i)  );
  //      }
}

void BlockMatrixModule::create_bases_functions(double bc_pc,const VecDouble &previous_sw,const VecDouble &previous_si,double por,double k,VecDouble  &newSi,VecDouble &newHomSi) 
{
//Initialize the vectors to contain the Pc from the fractures
// _SwIntegral.reinit(n_fracture_cells);
  
  assembleSystem(previous_sw, previous_si,k, por, bc_pc, _dt);
  // store newSi and newHomSi in m_ MSi and m_HomMSi
   newSi=Sol;
   newHomSi=Sol_hom;
}

double BlockMatrixModule::get_source(unsigned blockId,double new_bc_pc) 
{
  VecDoubleRef PSi(&(m_MSi(blockId,0)),_m);
  VecDoubleRef PSiHom(&(m_MHomSi(blockId,0)),_m);
  
  _BcPcNew(blockId)=new_bc_pc;
  
  VecDouble Sol;

  //Sol = PSi + (PcNew - PcOld) *PSiHom
  Sol=PSi;
  
  Sol.sadd(1,new_bc_pc - _BcPcOld(blockId),PSiHom);

  //Sol =  PcInv(Sol)
  _fpcinv.evaluate(Sol);
  
  //Now Sol vector contains  the saturation
  double SwInt=_mesh.getIntegralAtVertices(Sol);
  return -(SwInt-_SwIntegral(blockId))*_vPor(blockId)/(_dt*_Bx);
}

BlockMatrixModule::BlockMatrixModule(OrthoMesh &mesh,unsigned n_fracture_cell,  VecDouble &vPcInitial,VecDouble &vPor, VecDouble &vK, Function1D &fdpc, Function1D &fpcinv, Function1D &fmobt, Function1D &fmobrW): _mesh(mesh),n_fracture_cells(n_fracture_cell),_BcPcOld(vPcInitial), _vPor(vPor), _vK(vK) ,_fdpc(fdpc), _fpcinv(fpcinv), _fmobt(fmobt), _fmobrW(fmobrW)
{
  _dt=1.0;
  _dx = _mesh.get_dx();
  _dy = _mesh.get_dy();
  _dz = _mesh.get_dz();
  //  vPcInitial.print();
  // cout<<"dx \t"<<_dx<<_dy<<_dz<<endl;
  _m =_mesh.numRawVertices();
  //Initialize internal structures
  m_MSi.reinit(n_fracture_cells,_m);
  m_MHomSi.reinit(n_fracture_cells,_m);
  //Initialize the vectors to contain the Pc from the fractures
  _SwIntegral.reinit(n_fracture_cells);
  _BcPcNew=vPcInitial; // Initial Pc

  _Bx = _mesh.meshVolume(); // volume of the matrix block   
  m_MHomSi=0;
  for (unsigned i=0; i < n_fracture_cells; i++)
  {
    m_MSi.fill_row(i,_BcPcOld(i));
  }
  sp.reinit(2*_m, 2*_m,14);
  rhs_Si.reinit(2*_m);
  rhs_Sihom.reinit(2*_m);
  Sol.reinit(_m);
  Sol_hom.reinit(_m);
  Sol_old.reinit(_m);
  Sol_hom_old.reinit(_m);
  sat_old.reinit(_m);
  Solution1.reinit(2*_m);
  Solution2.reinit(2*_m);
  buildSparsityPattern();

}

---