ktmethoddbyd.cpp

Go to the documentation of this file.

#include "ktmethoddbyd.h"





KTMethodDByD::KTMethodDByD(OrthoMesh &mesh,Function3D &fInitU,const VecDouble &cPor,Function3D &fPrescribedU,FaceFluxFunction &flux, FixedValueCondition  &fixedC,double CFL) 
  :LaxFriedrichsForSystem(mesh,fInitU,cPor,fPrescribedU,flux,fixedC,CFL)
 {
   DX_2[0]=mesh.get_dx()/2.0;
   DX_2[1]=mesh.get_dy()/2.0;
   DX_2[2]=mesh.get_dz()/2.0;
   assert(n_components == 1);
   m_cValues.setRef(&(m_sol(0,0)),m_sol.size());
   m_cValuesPrev.setRef(&(m_solAnt(0,0)),m_sol.size());

   
   _MG.reinit(mesh.numCells(),3);
   m_cValuesAux.reinit(mesh.numCells());
   

   return;

   
}


void KTMethodDByD::iterate(double t, double tEnd) 
{

  double timeInterval = tEnd-t;
  double dt = this->calculateTimeStep(m_mesh,timeInterval,m_CFL,getPorosity(),m_flux);
  int nIteraction = (int) round(timeInterval/dt); //number of iteractions

  printf("KT dt = %g\n",dt);
  //For each time step
  for (int count = 0; count < nIteraction; count++)
  {
    iterateOnce(m_cValuesPrev,m_cValues,dt);
    iterateOnce(m_cValuesAux,m_cValues,dt);
    m_cValues.sadd(0.5,0.5,m_cValuesPrev);
  }
}


void KTMethodDByD::iterateOnce(VecDouble &vAux, VecDouble &vSol, double dt) 
{
  buildDerivatives(m_mesh,vSol,_MG);


    //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();
  OrthoMesh::Face_It face = m_mesh.begin_face();
  OrthoMesh::Face_It endf = m_mesh.end_face();
  m_vbc.rewind();
  bool hasbc;

  vAux=vSol;

    //For each face
    for (;face!=endf;face++)
    {
      //Get the values of the previous solution in the right and left of the face
      static VecDouble SwNeg(1),SwPos(1),vFlux(1);
      unsigned negCell,posCell;

      face->getAdjCellIndices(negCell,posCell);

      //Get the values of the previous solution in the right and left of the face
      hasbc=this->getValuesOfTheFaceCellsScalar(m_mesh,m_vbc,vAux,face->index(),negCell,posCell,&SwNeg(0),&SwPos(0));


      //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;



      unsigned normal = face->getNormalOrientation();
      
      double posPor = getPorosity()(posCell);
      double negPor = getPorosity()(negCell);
      //double mPor = (posPor + negPor)/2.0;
      
      SwNeg(0) += DX_2[normal]*_MG(negCell,normal); 
      SwPos(0) -= DX_2[normal]*_MG(posCell,normal); 

      //calculate the right coefficient of diffusion 
      double a;
      a = fabs(m_flux.maxLocalCharVelocity(face,SwNeg,SwPos));

      //printf("%g ",a);
       //Calculate the amount of the mass passed through the face and divide it by the cell volume
      m_flux.fluxAtFace(vFlux,face,SwNeg,SwPos);
      double flux =  dt*face->areaPerCellVol()* (vFlux(0)  - a*(SwPos(0)-SwNeg(0))/2.0);
      //printf("face %d = %g,Pos: %d,Neg %d, SwPos = %g,SwNeg = %g \n",faceIndex,flux,face->getPosCellIndex(),face->getNegCellIndex(),SwPos,SwNeg);
      if (face->hasNegCell())
        vSol(negCell) +=  - flux/negPor;

      if (face->hasPosCell())
        vSol(posCell) +=  + flux/posPor;
    }
}

---