laxfriedrichsforsystem.cpp

Go to the documentation of this file.

#include "laxfriedrichsforsystem.h"
#include "numericmethods.h"
#include "hdf5orthowriter.h"
#include "gnuplotanim.h"
#include  "flashcompositional.h"
#include "dynamicbase.h"

LaxFriedrichsForSystem::LaxFriedrichsForSystem(OrthoMesh &mesh,Function3D &fInitU,const VecDouble &cPor,Function3D &fPrescribedU,FaceFluxFunction &flux, FixedValueCondition &fixedC,double CFL) 
  :m_mesh(mesh),m_cPor(cPor),m_flux(flux),n_components(flux.numComponents()),m_fixedC(fixedC)
    

 {
   assert(fInitU.n_components() == n_components);
   assert(fPrescribedU.n_components() == n_components);
  assert(cPor.size() == mesh.numCells());
  m_CFL = CFL;
  

  //Alloc Vectors
  m_sol.reinit(mesh.numCells(),n_components);
  m_solAnt.reinit(mesh.numCells(),n_components);

  
  
  mesh.setCentralValuesFromFunction(fInitU,m_sol);
    

  //Initialize the boundary condtions
  m_vbc.set_bc_from_function(mesh,fPrescribedU);

  //  std::cout << m_vbc << std::endl;
  //exit(0);
  //Initial conditions
  //m_vSol.reinit(mesh.numCells(),n_components);
  //m_vSolAnt.reinit(mesh.numCells(),n_components);
  //mesh.setCentralValuesFromFunction(fInitU,m_vSol);
  
 }



ArrayOfVecDouble& LaxFriedrichsForSystem::getSolutionAtCells() 
{
  return m_sol;
}



void LaxFriedrichsForSystem::updateVelocities(DynamicBase &dynMod) 
{
  m_flux.updateDynamicData(dynMod);
  FlashCompositional *flash = getFlash();

  if (flash)
    this->updateBCForCompressibility(m_vbc,dynMod,*flash);
  
}



double LaxFriedrichsForSystem::getDt(double t,double tEnd) 
{
  assert(tEnd >= t);
  double timeInterval = tEnd-t;
  double dt = this->calculateTimeStep(m_mesh,timeInterval,m_CFL,getPorosity(),m_flux);
  return dt;
}


void LaxFriedrichsForSystem::iterateN(unsigned nSteps, double dt) 
{
  unsigned posCell,negCell;

  //We iterate the method along faces, so 
  //alloc the vectors to contain the solutions in neg
  //and pos cells of the faces.  Remeber that the solution
  //in a cell is vector with a value for each component
  static VecDouble  vBC(n_components),vFlux(n_components);
  static VecDoubleRef vQPos(n_components),vQNeg(n_components);
  //For each time step
  for (unsigned count = 0; count < nSteps; count++)
  {
    //Get the face iterator and the number of faces.
    OrthoMesh::Cash_Face_It face = m_mesh.begin_cash_face();
    OrthoMesh::Cash_Face_It endf = m_mesh.end_cash_face();
    m_vbc.rewind();
    bool hasbc;

    m_solAnt=m_sol;
    
    //For each face
    for (;face!=endf;face++)
    {
      face->getAdjCellIndices(negCell,posCell);
      //Get the values of the previous solution in the right and left of the face
      hasbc=this->getValuesOfTheFaceCells(m_mesh,m_vbc,m_solAnt,face->index(),negCell,posCell,vQNeg,vQPos);

      //printf("face %d : Values %g %g %g | %g %g %g\n",face->index(),vQNeg(0),vQNeg(1),vQNeg(2),vQPos(0),vQPos(1),vQPos(2));
      
      //Get the indices of the cells that contain the face.
      //If the face is at boundary, it has just one cell.
      //In this case the method getAdjCells() return posCell == negCell
      if (posCell == OrthoMesh::INVALID_INDEX)
        posCell=negCell;
      else if (negCell == OrthoMesh::INVALID_INDEX)
        negCell=posCell;
        


      double posPor = getPorosity()(posCell);
      double negPor = getPorosity()(negCell);
      double mPor = (posPor + negPor)/2.0;


      //get Flux
      m_flux.fluxAtFace(vFlux,face,vQNeg,vQPos);
      
      //Now update the cell values component by component
      for (unsigned cmp=0;cmp<n_components;cmp++)
      {
        
        //See if the face has prescribed boundary condition BC      
        double flux;
        if (!hasbc)    
        {
          /*
            The face doesnt have boundary condition BC
            proceed to the usual computation of the flux
            to calculate the amount of the mass passed through
            the face and divide it by the cell volume
          */
           
          flux = dt*face->areaPerCellVol() * vFlux(cmp) - mPor*(vQPos(cmp)-vQNeg(cmp))/2.0;
        }
        else //the face has BC, so dont add the stability term
        {
          flux = dt*face->areaPerCellVol() * vFlux(cmp);
        }

      
      
        //printf("face %d = %g,Pos: %d,Neg %d, SwPos = %g,SwNeg = %g \n",faceIndex,flux,face->getPosCellIndex(),face->getNegCellIndex(),SwPos,SwNeg);
        if (face->hasPosCell())
          m_sol(posCell,cmp) +=  + flux/posPor;

        if (face->hasNegCell())
          m_sol(negCell,cmp) +=  - flux/negPor;
        
      } //end for each component
    } //end for each face
  } //end for each time step
}






 void LaxFriedrichsForSystem::printOutput() 
{
  VecDouble vValues(m_mesh.numVertices());
  VecDouble cValues(m_mesh.numCells());

  double acum=0.0;
  for (unsigned cmp=0;cmp<n_components;cmp++)
  {
    char str[100];
    m_sol.getComponent(cValues,cmp);
    this->m_mesh.projectCentralValuesAtVertices(cValues,vValues);
    //printf("Mass S%d: %g\n",cmp+1,m_mesh.getIntegralAtCells(cValues);
    acum+=m_mesh.getIntegralAtCells(cValues);
    HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 
    sprintf(str,"S%d",cmp+1);
    hdf5.writeScalarField(vValues,str);

  }
}




 void LaxFriedrichsForSystem::updateFixedPressureBC(Function3D *fDP,FlashCompositional &flash) 
{
  FlashData flashData(flash.numPhases(),flash.numComponents());
  VecDouble phasesVol(flash.numPhases());
  assert(n_components == flash.numComponents());
  OrthoMesh::Face_It fend = m_mesh.end_face();
  flashData.allocateOwnMemory();
  Point3D pt;
  VecDoubleRef bc(n_components);
  m_vbc.rewind();
  //for all faces
  //m_fBC[0].print(std::cout);
  for (OrthoMesh::Face_It face=m_mesh.begin_face(); face!=fend; face++)
    {
      //If face is at boundary and there is boundary conditions for the transport
      unsigned i = face->index();
      if (m_vbc.get_ref_bc(i,bc))
      {
        face->barycenter(pt);
        if (fDP->isInDomain(pt))
        {
          //Get components for face i

          double P = (*fDP)(pt);

          //Make the flash and get the phase volumes
          flash.flash(P,bc,flashData);
          flash.getPhasesVolume(P,flashData,phasesVol);

          //calculate factor
          double factor=1.0/NumericMethods::vectorSum(phasesVol);
          //printf("%g %g %g\n",factor,phasesVol(0),phasesVol(1));
          //flashData.print();
          bc*=factor;
        }
      }
    }
  std::cout << m_vbc;
}









void LaxFriedrichsForSystem::updateBCForCompressibility(VecIncBC &vbc,DynamicBase &dyn,FlashCompositional &flash) 
{
  const VecDouble &P = dyn.getPressureAtCells();
  VecIncBC::data_iterator itBC = vbc.begin_data();
  VecDouble phasesVol(flash.numPhases());
  FlashData flashData(flash.numPhases(),flash.numComponents());
  VecDoubleRef m(vbc.bc_dim());
  for (;itBC!= vbc.end_data();itBC++)
  {
    //get the bc
    m.setData(itBC->data.p);
    //appy the flash
    flash.flash(P(itBC->data.nearCell),m,flashData);
    //calculate the volumes
    flash.getPhasesVolume(P(itBC->data.nearCell),flashData,phasesVol);
    double factor=1.0/NumericMethods::vectorSum(phasesVol);
    //Multiply the components by such factor
    m*=factor;
  }
}

---