upwindforsystem.cpp

Go to the documentation of this file.

#include "upwindforsystem.h"
#include "numericmethods.h"


 UpwindForSystem::~UpwindForSystem() 
{

}

UpwindForSystem::UpwindForSystem(OrthoMesh &mesh,Function3D &fInitU,const VecDouble &cPor,Function3D &fPrescribedU,FaceFluxFunction &flux, FixedValueCondition &fixedC,double CFL,DiffusiveStep *diffStp) 
  :LaxFriedrichsForSystem(mesh,fInitU,cPor,fPrescribedU,flux,fixedC,CFL)
{

}




void UpwindForSystem::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
  VecDouble  vBC(n_components),vFluxPos(n_components),vFluxNeg(n_components);
  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::Face_It face = m_mesh.begin_face();
    OrthoMesh::Face_It endf = m_mesh.end_face();


    //m_solAnt receives the solution in time tn since
    //m_sol will contain the solution in tn+1. 
    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
      this->getValuesOfTheFaceCells(m_mesh,m_vbc,m_solAnt,face->index(),negCell,posCell,vQNeg,vQPos);
      
      //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;

      //Get the porsity
      double posPor = getPorosity()(posCell);
      double negPor = getPorosity()(negCell);


      m_flux.fluxAtFace(vFluxPos,face,vQPos,vQPos);
      m_flux.fluxAtFace(vFluxNeg,face,vQNeg,vQNeg);

      //Now get the flux component by component
      for (unsigned cmp=0;cmp<n_components;cmp++)
      {
        assert(vQPos(cmp) >= 0.0 && vQNeg(cmp) >= 0.0);

        double fluxPos,fluxNeg,flux;

        //Get the flux from the positive and negative side
        fluxPos = dt*face->areaPerCellVol() * vFluxPos(cmp);
        fluxNeg = dt*face->areaPerCellVol() * vFluxNeg(cmp);

        
        //if the flux is coming from the positive cell, use the flux computed from the value of the  Pos Cell
        if (fluxPos >= 0.0 && fluxNeg >= 0.0)
          flux = fluxNeg;
        else if (fluxPos <0.0 && fluxNeg < 0.0) //if the flux comes from the other direction use the Neg flux
          flux = fluxPos;
        else
        {
          flux = (fluxPos + fluxNeg)/2.0 - (posPor+negPor)/2.0 * (vQPos(cmp)-vQNeg(cmp))/6.0; 
          //printf("Sonic Point: %g %g %d\n",fluxPos,fluxNeg,face->getNormalOrientation());
        }
        if (OrthoMesh::isValid(posCell))
          m_sol(posCell,cmp) +=  + flux/posPor;

        if (OrthoMesh::isValid(negCell))
          m_sol(posCell,cmp) +=  - flux/negPor;
        
      } //end for each component
    } //end for each face
  } //end for each time step

}

---