compdiffusivestep.cpp

Go to the documentation of this file.

#include "compdiffusivestep.h"
#include "numericmethods.h"
#include "hdf5orthowriter.h"

using NumericMethods::_div;

CompDiffusiveStep::CompDiffusiveStep(OrthoMesh &mesh, Function3D &fPrescribedPressure, Function3D &fTransportBC,const VecDouble &K, const VecDouble &vPor,Function1D &fMobPc, Function1D &DfPc,FlashCompositional &flash,LinearSolver &solver) 
  :m_mesh(mesh),m_solver(solver),m_flash(flash),m_vK(K),m_vPor(vPor),m_fMobPc(fMobPc),m_DfPc(DfPc)
 {
   m_vCoeff.reinit(mesh.numCells());
   m_vCoeff2.reinit(mesh.numCells());
   m_sol.reinit(mesh.numCells());
   m_B.reinit(mesh.numCells());

   m_vSaq.reinit(mesh.numCells());
   m_vSaqPrev.reinit(mesh.numCells());
   
   //setSaqBC(m_flashBC,m_swBC,fPrescribedPressure,fTransportBC,flash);
   
   buildSparsityPattern(mesh,m_sparsePattern);
   m_A.reinit(m_sparsePattern);
   m_cV.reinit(mesh.numCells());
   m_vV.reinit(mesh.numVertices());

 }

void CompDiffusiveStep::setSaqBC(VecIncAccess<FlashBC> &flashbc, ScalarBC &swbc,const Function3D &fPressureBC, const Function3D &fTransBC,FlashCompositional &flash) 
{
  assert(flash.numPhases() == 2);
  assert(fTransBC.getImageDim() == flash.numComponents());
  unsigned ncmps = flash.numComponents();
  
  VecDouble cmps(flash.numComponents());
  VecDouble phasesVol(flash.numPhases());
  FlashData flashData(flash.numPhases(),flash.numComponents());
  flashData.allocateOwnMemory();

  OrthoMesh::Face_It face = m_mesh.begin_face();
  OrthoMesh::Face_It endf = m_mesh.end_face();
  for (;face!=endf;face++)
  {
    if (face->at_boundary())
    {
      Point3D p = face->barycenter();
      if (fTransBC.isInDomain(p))
      {
        assert(fPressureBC.isInDomain(p));
        for (unsigned i=0;i<ncmps;i++)
        {
          cmps(i)=fTransBC(p,i);
        }
        double P = fPressureBC(p);
        flash.flash(P,cmps,flashData);
        flash.getPhasesVolume(P,flashData,phasesVol);
        double Saq = phasesVol(0)/(phasesVol(1)+phasesVol(0));
        FlashBC  _flashBC; //(flash.numPhases(),flash.numComponents());;
        _ScalarBC _sBC;
        _flashBC.Saq=Saq;
        _sBC.value=Saq;
        _sBC.face=face;
        swbc.push_back(_sBC);
        flashbc.push_back(face->index(),_flashBC);
      }
    }
  }
}



void CompDiffusiveStep::iterate(double dt) 
{
  const VecDouble &vP=m_flash.getDynamicModule().getPressureAtCells();
  TransportBase &trans = m_flash.getTransportModule();
  assert(m_flash.numPhases()==2);
  
  static VecDouble phasesVol(m_flash.numPhases());




  //Allocate the flash data
  static FlashData data(m_flash.numPhases(),m_flash.numComponents());


  Index nCells = m_mesh.numCells();
  double cellVol=m_mesh.begin_cell()->volume();  

  //Write Zeros in the system
  m_A=0.0;
  m_B=0.0;

  for (Index cellIndex=0;cellIndex < nCells; cellIndex++)
  {
    //Compute/Get the flash data
    m_flash.flash(cellIndex,data);
    double P=vP(cellIndex);
    m_flash.getPhasesVolume(P,data,phasesVol);     


    double Saq=phasesVol(0)/(phasesVol(1)+phasesVol(0));
    m_vSaqPrev(cellIndex)=Saq;
    double por=m_vPor(cellIndex);
    //double molarVol = Saq/data.getPhaseTotalMoles(0);    
    double molarVol;
    if (fabs(Saq) == 0.0)
      molarVol=1.2733e-4;
    else
      molarVol = Saq/data.getPhaseTotalMoles(0);
    
    m_A.add(cellIndex,cellIndex,molarVol*por*cellVol/dt);
    m_B(cellIndex)+=molarVol*por*Saq/dt*cellVol;

    m_vCoeff(cellIndex)=-m_vK(cellIndex)*m_DfPc(Saq)*m_fMobPc(Saq);
    m_vCoeff2(cellIndex)=m_vCoeff(cellIndex)*molarVol;
  }

  //Now solve the system
  MixedHybridBase::assemblySystem(m_A,m_B,m_mesh,m_swBC,m_vCoeff2);


  m_solver.solve(m_A,m_vSaq,m_B);
  unsigned n_components = m_flash.numComponents();

  ArrayOfVecDouble &vSol = trans.getSolutionAtCells();

  //m_sol contains the saturation


  for (unsigned cmp=0;cmp<n_components;cmp++)
  {
    OrthoMesh::Face_It face = m_mesh.begin_face();
    OrthoMesh::Face_It end = m_mesh.end_face();
    Index c1,c2;
    FlashData data1(m_flash.numPhases(),m_flash.numComponents());
    FlashData data2(m_flash.numPhases(),m_flash.numComponents());
    VecIncAccess<FlashBC>::iterator it=m_flashBC.begin();
    for (;face!=end;face++)
    {
      if (!face->at_boundary())
      {
        face->getAdjCellIndices(c1,c2);
        double Saq1=m_vSaqPrev(c1);
        double Saq2=m_vSaqPrev(c2);
        m_flash.flash(c1,data1);
        m_flash.flash(c2,data2);
        double Kflux1 = m_vCoeff(c1)*(_div(data1.getMoles(0,cmp),Saq1) -  _div(data1.getMoles(1,cmp),(1-Saq1)));
        double Kflux2 = m_vCoeff(c2)*( _div(data2.getMoles(0,cmp),Saq2) - _div(data2.getMoles(1,cmp),(1-Saq2)));
        double flux =  - 2*_div(Kflux1*Kflux2, (Kflux1+Kflux2))*(m_vSaq(c2)-m_vSaq(c1))*face->areaPerCellVol();
        //double flux =  -(Kflux1+Kflux2)/2.0*(m_vSaq(c2)-m_vSaq(c1))*face->areaPerCellVol();

        //if (cmp == 0)
        //  printf("%d) %g %g %g\n",face->index(),flux,Kflux1,Kflux2);
        vSol(c1,cmp)-=flux/m_vPor(c1)*dt*face->areaPerCellVol();
        vSol(c2,cmp)+=flux/m_vPor(c2)*dt*face->areaPerCellVol();
      }
      else
      {
        Index faceIndex = face->index();
        const FlashBC *flash_bc = it.inc_search(faceIndex);
        if (flash_bc)
        {
          abort();
          face->getAdjCellIndices(c1,c2);
          double Saq = flash_bc->Saq;
          if (c1 != OrthoMesh::INVALID_INDEX)
          {  
            m_flash.flash(c1,data1);
            double Kflux = m_vCoeff(c1)*(_div(data1.getMoles(0,cmp),m_vSaqPrev(c1)) - _div(data1.getMoles(1,cmp),1-m_vSaqPrev(c1)));
            double flux= -2*Kflux*(Saq-m_vSaq(c1))*face->areaPerCellVol();
            vSol(c1,cmp)-=flux/m_vPor(c1)*dt*face->areaPerCellVol();
          }
          else
          {
            m_flash.flash(c2,data2);
            double Kflux = m_vCoeff(c2)*(_div(data2.getMoles(0,cmp),m_vSaqPrev(c2)) - _div(data2.getMoles(1,cmp),1-m_vSaqPrev(c2)));
            double flux= -2*Kflux*(m_vSaq(c2) - Saq)*face->areaPerCellVol();
            vSol(c2,cmp)+=flux/m_vPor(c2)*dt*face->areaPerCellVol();
          }
        }
      }
    }
  }
}


 CompDiffusiveStep::~CompDiffusiveStep() 
{

}

---