conservativemethodforsystem.cpp

Go to the documentation of this file.

#include "conservativemethodforsystem.h"
#include "numericmethods.h"
#include "flashbase.h"
#include "vecincbc.h"
#include "arrayofvecdouble.h"

double ConservativeMethodForSystem::calculateTimeStep(OrthoMesh &mesh,double tEnd,double MaxCFL, const VecDouble &vPorosity, FaceFluxFunction &flux) 
{
  double dMinPor = NumericMethods::vectorMinValue(vPorosity);
  double dInvMinTravelTime; //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);

  //To avoid division by zero we calculate the inverse
  //of the min travel time and then we invert the number.

  //Min time to cross a cell  in the X direction
  dInvMinTravelTime = fabs(MaxCharVel[X])/dX;

  //Min time to cross a cell  in the Y direction
  dInvMinTravelTime = fmax(fabs(MaxCharVel[Y])/dY, dInvMinTravelTime);

 //Min time to cross a cell  in the Z direction
 dInvMinTravelTime = fmax(fabs(MaxCharVel[Z])/dZ, dInvMinTravelTime);

 
 if (dInvMinTravelTime == 0)
   return tEnd;

 //double dTravelTime = 1.0/dInvMinTravelTime;
 
 if (tEnd == INFINITY)
   return dMinPor*MaxCFL/dInvMinTravelTime; 
 else
   return tEnd/ceil(tEnd*dInvMinTravelTime/(dMinPor*MaxCFL));

}









bool ConservativeMethodForSystem::getValuesOfTheFaceCells(OrthoMesh &mesh,VecIncBC &vbc,const ArrayOfVecDouble& cvcValues, Index face_index, Index negCell, Index posCell,VecDoubleRef &vSwNeg,VecDoubleRef &vSwPos) 
{
  ArrayOfVecDouble &vcValues=const_cast<ArrayOfVecDouble&>(cvcValues);
  assert(vSwNeg.size() == vbc.bc_dim());
  assert(vSwPos.size() == vbc.bc_dim());
  if (posCell != OrthoMesh::invalidIndex())
  {
    //get the positive value
    vcValues.getVecValues(posCell,&vSwPos);

    //get the negative value if it exist
    if (negCell != OrthoMesh::invalidIndex())
    {
      vcValues.getVecValues(negCell,&vSwNeg);
      return false;
    }
    else
    {

      //if the neg value doesnt exist,
      //try to see if there is boundary condition
      //for the face
      if (vbc.get_ref_bc(face_index,vSwNeg))
      {
        return true;
      }
      else
      {
        //If not use mirror condition
        vcValues.getVecValues(posCell,&vSwNeg);
        return false;
      }
   }
  }
  else
  {
    //The positive value does not exist,
    //So the Negative value must exist. Get it
    vcValues.getVecValues(negCell,&vSwNeg);

    //Now take the value for SwPos from the boundary
    //condition or make it equal to the value of the negative cell.
    if (vbc.get_ref_bc(face_index,vSwPos)) //Ok BC exists
    {
      return true;
    }
    else
    {
      //If not use mirror condition
      vcValues.getVecValues(negCell,&vSwPos);
      return false;
    }
  }
  assert(0);
  return false;
}

bool ConservativeMethodForSystem::getValuesOfTheFaceCellsScalar(OrthoMesh &mesh,VecIncBC &vbc,const VecDouble &cValues, Index face_index, Index negCell, Index posCell,double *dNeg,double *dPos) 
{
  bool result=false;
  static VecDouble vSw(1);
  assert(vbc.bc_dim()==1);
  if (posCell != OrthoMesh::invalidIndex())
  {
    //get the positive value
    *dPos = cValues(posCell);

    //get the negative value if it exist
    if (negCell != OrthoMesh::invalidIndex())
      *dNeg=cValues(negCell);
    else
    {

      //if the neg value doesnt exist,
      //try to see if there is boundary condition
      //for the face
      if (vbc.copy_bc(face_index,vSw))
      {
        result=true;
        *dNeg=vSw(0);
      }
      else
      {
        //If not use mirror condition
        *dNeg=*dPos;
      }
   }
  }
  else
  {
    //The positive value does not exist,
    //So the Negative value must exist. Get it
    *dNeg = cValues(negCell);

    //Now take the value for SwPos from the boundary
    //condition or make it equal to the value of the negative cell.
    if (vbc.copy_bc(face_index,vSw))
    {
      *dPos=vSw(0);
      result=true;
    }
    else
    {
      *dPos=*dNeg; 
    }
  }
  return result;
  
}




void ConservativeMethodForSystem::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 ConservativeMethodForSystem::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();
}

---