conservativemethod.cpp

Go to the documentation of this file.

#include "conservativemethod.h"
#include "numericmethods.h"


void ConservativeMethod::getValuesOfTheFaceCells(OrthoMesh &mesh,const VecDouble &fBCValues,const VecDouble& cValues,Index face_index, Index cNeg,Index cPos,double &SwNeg,double &SwPos)
{
  assert(fBCValues.size() == mesh.numFaces());
  assert(cValues.size() == mesh.numCells());

  static unsigned Null = OrthoMesh::invalidIndex();
  if (cNeg != Null)
  {
    //get the negative value
    SwNeg = cValues(cNeg);

    //get the positive value if it exist
    if (cPos != Null)
      SwPos = cValues(cPos);
    else
    {

      //if the pos value doesnt exist,
      //try to see if there is boundary condition
      //for the face
      SwPos = fBCValues(face_index);
      if (SwPos == INFINITY) //boundary condition is not defined ?
        SwPos=SwNeg; //If not exist boundary condition, use the value of the negative cell
    }
  }
  else
  {
    //The negative  value does not exist,
    //So the Positive value must exist. Get it
    SwPos = cValues(cPos);

    //Now take the value for SwNeg from the boundary
    //condition or make it equal to the value of the positive cell.
    SwNeg = fBCValues(face_index);
    if (SwNeg == INFINITY) //Is boundary condition not defined ?
      SwNeg = SwPos;
  }
}


void ConservativeMethod::getValuesOfTheFaceCells(OrthoMesh &mesh,const VecDouble &fBCValues,const VecDouble& cValues,OrthoMesh::Cash_Face_It & face,double &SwNeg,double &SwPos)
{
  assert(fBCValues.size() == mesh.numFaces());
  assert(cValues.size() == mesh.numCells());
  unsigned cPos,cNeg;
  face->getAdjCellIndices(cNeg,cPos);

  unsigned Null = OrthoMesh::invalidIndex();
  if (cPos != Null)
  {
    //get the positive value
    SwPos = cValues(cPos);
    //get the negative value if it exist
    if (cNeg != Null)
      SwNeg = cValues(cNeg);
    else
    {

      //if the neg value doesnt exist,
      //try to see if there is boundary condition
      //for the face
      SwNeg = fBCValues(face->index());
      if (SwNeg == INFINITY) //boundary condition is not defined ?
        SwNeg = SwPos; //If not exist boundary condition, use the value of the positive cell
    }
  }
  else
  {
    //The positive value does not exist,
    //So the Negative value must exist. Get it
    SwNeg = cValues(cNeg);

    //Now take the value for SwPos from the boundary
    //condition or make it equal to the value of the negative cell.
    SwPos = fBCValues(face->index());
    if (SwPos == INFINITY) //Is boundary condition not defined ?
      SwPos = SwNeg;
  }
}


void ConservativeMethod::getValuesOfTheFaceCells(OrthoMesh &mesh,const VecDouble &fBCValues,const VecDouble& cValues,OrthoMesh::Face_It & face,double &SwNeg,double &SwPos)
{
  assert(fBCValues.size() == mesh.numFaces());
  assert(cValues.size() == mesh.numCells());
  unsigned cPos,cNeg;
  face->getAdjCellIndices(cNeg,cPos);

  unsigned Null = OrthoMesh::invalidIndex();
  if (cPos != Null)
  {
    //get the positive value
    SwPos = cValues(cPos);
    //get the negative value if it exist
    if (cNeg != Null)
      SwNeg = cValues(cNeg);
    else
    {

      //if the neg value doesnt exist,
      //try to see if there is boundary condition
      //for the face
      SwNeg = fBCValues(face->index());
      if (SwNeg == INFINITY) //boundary condition is not defined ?
        SwNeg = SwPos; //If not exist boundary condition, use the value of the positive cell
    }
  }
  else
  {
    //The positive value does not exist,
    //So the Negative value must exist. Get it
    SwNeg = cValues(cNeg);

    //Now take the value for SwPos from the boundary
    //condition or make it equal to the value of the negative cell.
    SwPos = fBCValues(face->index());
    if (SwPos == INFINITY) //Is boundary condition not defined ?
      SwPos = SwNeg;
  }
}




double  ConservativeMethod::calculateTimeStep(OrthoMesh &mesh,double tEnd,double MaxCFL, const VecDouble &vPorosity, FaceFluxFunction &flux)
{
  double dMinPor = NumericMethods::vectorMinValue(vPorosity);
  double dMinTravelTime; //Travel time is the time the solution needs to cross a cell.

  //Get the discretization steps in each dimension
  double dX = mesh.getDX()[0];
  double dY = mesh.getDX()[1];
  double dZ = mesh.getDX()[2];

  //Get the maximum chracteristic velocity in the three components
  double MaxCharVel[3];
  flux.maxGlobalCharVelocity(MaxCharVel);
  dMinTravelTime=0.0;

  //Min time to cross a cell  in the X direction
  if (MaxCharVel[X] != 0.0)
    dMinTravelTime = dX/fabs(MaxCharVel[X]);
  printf("Max Vel %g %g %g\n",MaxCharVel[0],MaxCharVel[1],MaxCharVel[2]);

  //Min time to cross a cell  in the Y direction
  if (MaxCharVel[Y] != 0.0)
    dMinTravelTime = fmin(dY/fabs(MaxCharVel[Y]), dMinTravelTime);

 //Min time to cross a cell  in the Z direction
  if (MaxCharVel[Z] != 0.0)
    dMinTravelTime = fmin(dZ/fabs(MaxCharVel[Z]), dMinTravelTime);

 return tEnd/ceil(tEnd/(dMinPor*dMinTravelTime*MaxCFL));
 
}



void ConservativeMethod::buildDerivatives(OrthoMesh &mesh,const VecDouble &cV,Matrix &MG) 
{
 
  assert(MG.m() == mesh.numCells());
  assert(MG.n() == 3);
  assert(cV.size() == mesh.numCells());

  OrthoMesh::Cell_It cell = mesh.begin_cell();
  OrthoMesh::Cell_It endc = mesh.end_cell();
  double s0,s1,s2;
  for (;cell != endc; cell++)
  {
    getCellDerivatives(cV,cell,&s0,&s1,&s2);
    MG(cell->index(),0)=s0;
    MG(cell->index(),1)=s1;
    MG(cell->index(),2)=s2;
  } 
}



void ConservativeMethod::getCellDerivatives(const VecDouble &sol,const OrthoMesh::Cell_It &cell,double *s0,double *s1, double *s2) 
{
  //Derivative for x
  const OrthoMesh &mesh=cell->getMesh();
  Index li,ri;
  double u = sol(cell->index());
  li = cell->neighbor_index(LEFT_CELL);
  ri = cell->neighbor_index(RIGHT_CELL);
  if (li == OrthoMesh::INVALID_INDEX || ri == OrthoMesh::INVALID_INDEX)
     *s0 = 0.0;
  else
    *s0=NumericMethods::minMod(u-sol(li),sol(ri)-u)/mesh.get_dx();
  

  li = cell->neighbor_index(FRONT_CELL);
  ri = cell->neighbor_index(BACK_CELL);
  if (li == OrthoMesh::INVALID_INDEX || ri == OrthoMesh::INVALID_INDEX)
     *s1 = 0.0;
  else
    *s1=NumericMethods::minMod(u-sol(li),sol(ri)-u)/mesh.get_dy();

  li = cell->neighbor_index(BOTTOM_CELL);
  ri = cell->neighbor_index(UP_CELL);
  if (li == OrthoMesh::INVALID_INDEX || ri == OrthoMesh::INVALID_INDEX)
     *s2 = 0.0;
  else
    *s2=NumericMethods::minMod(u-sol(li),sol(ri)-u)/mesh.get_dz();
}

---