laxfriedrichssystemmpi.cpp

Go to the documentation of this file.

#include "laxfriedrichssystemmpi.h"
#include "netmpi.h"
#include "fxyslabid.h"
#define SEND_GHOSTCELL 201
#define DEAD_ZONE_VALUE 202



LaxFriedrichsSystemMPI::LaxFriedrichsSystemMPI(OrthoMesh &mesh,Function3D &fInitU,const VecDouble &cPor,Function3D &fPrescribedU,FaceFluxFunction &flux, double CFL,int mesh_overlap) 
  :LaxFriedrichsForSystem(mesh,fInitU,cPor,fPrescribedU,flux,CFL)
{
  cLDeadZoneRcv.resize(n_components);
  cLDeadZoneSnd.resize(n_components);
  LGhostCellRcv.resize(n_components);
  LGhostCellSnd.resize(n_components);
  cRDeadZoneRcv.resize(n_components);
  cRDeadZoneSnd.resize(n_components);
  RGhostCellRcv.resize(n_components);
  RGhostCellSnd.resize(n_components);

  NetMPI::trace("System:: set Communication\n");
  //  if (NetMPI::rank() != 0)
  //  NetMPI::debugPoint();
  LaxFriedrichsSystemMPI::setCommunicationBuffers(mesh,mesh_overlap);
  NetMPI::trace("System:: end set Communication\n");
   
}


 LaxFriedrichsSystemMPI::~LaxFriedrichsSystemMPI() 
{

}


void LaxFriedrichsSystemMPI::setCommunicationBuffers(OrthoMesh &mesh,int mesh_overlap) 
{
  //Set communication buffers
  unsigned elemsInXY = mesh.numElemY()*mesh.numElemZ();
  double dx = mesh.get_dx();
  double tol = dx/4.0;

  fDeadZone.resize(mesh.numFaces(),ACTIVE_ZONE);

  if (NetMPI::rank() != 0)
  {
    double X0 = mesh.getP()[0];
    double X1 = X0 + mesh_overlap*dx;
    double X2 = X0 + 2*mesh_overlap*dx;
    
    cLDeadZoneRcvMap.reserve(mesh_overlap*elemsInXY);
    cLDeadZoneSndMap.reserve(mesh_overlap*elemsInXY);
    FXYSlabId fRcv(X0,X1,tol);
    FXYSlabId fSnd(X1,X2,tol);
    mesh.getCellsIndicesInDomain(fRcv,cLDeadZoneRcvMap);
    mesh.getCellsIndicesInDomain(fSnd,cLDeadZoneSndMap);


    LGhostCellRcvMap.reserve(elemsInXY);
    LGhostCellSndMap.reserve(elemsInXY);
    FXYSlabId f2Rcv(X1-dx,X1,tol);
    FXYSlabId f2Snd(X1,X1+dx,tol);
    mesh.getCellsIndicesInDomain(f2Rcv,LGhostCellRcvMap);
    mesh.getCellsIndicesInDomain(f2Snd,LGhostCellSndMap);
    NetMPI::trace("Start");

    assert(LGhostCellRcvMap.size() == elemsInXY);
    assert(LGhostCellSndMap.size() == elemsInXY);



    for (unsigned i=0;i<n_components;i++)
    {
      cLDeadZoneRcv[i].reinit(cLDeadZoneRcvMap.size());
      cLDeadZoneSnd[i].reinit(cLDeadZoneSndMap.size());
      LGhostCellRcv[i].reinit(LGhostCellRcvMap.size());
      LGhostCellSnd[i].reinit(LGhostCellSndMap.size());
    }

    FXYSlabId fLFaces(X0,X1-dx,tol);
    mesh.tagFacesInDomain(fLFaces,fDeadZone, LEFT_DEAD_ZONE);


  }
  if (!NetMPI::isLastProcess())
  {
    double X2 = mesh.getQ()[0];
    double X1 = X2 - mesh_overlap*dx;
    double X0 = X2 -  2*mesh_overlap*dx;

    cRDeadZoneRcvMap.reserve(mesh_overlap*elemsInXY);
    cRDeadZoneSndMap.reserve(mesh_overlap*elemsInXY);
    FXYSlabId fRcv(X1,X2,tol);
    FXYSlabId fSnd(X0,X1,tol);
    mesh.getCellsIndicesInDomain(fRcv,cRDeadZoneRcvMap);
    mesh.getCellsIndicesInDomain(fSnd,cRDeadZoneSndMap);

    RGhostCellRcvMap.reserve(elemsInXY);
    RGhostCellSndMap.reserve(elemsInXY);
    FXYSlabId f2Rcv(X1,X1+dx,tol);
    FXYSlabId f2Snd(X1-dx,X1,tol);
    mesh.getCellsIndicesInDomain(f2Rcv,RGhostCellRcvMap);
    mesh.getCellsIndicesInDomain(f2Snd,RGhostCellSndMap);
    NetMPI::trace("Start");
    assert(RGhostCellRcvMap.size() == elemsInXY);
    assert(RGhostCellSndMap.size() == elemsInXY);

    
    for (unsigned i=0;i<n_components;i++)
    {
      cRDeadZoneRcv[i].reinit(cRDeadZoneRcvMap.size());
      cRDeadZoneSnd[i].reinit(cRDeadZoneSndMap.size());
      RGhostCellRcv[i].reinit(RGhostCellRcvMap.size());
      RGhostCellSnd[i].reinit(RGhostCellSndMap.size());
    }

    FXYSlabId fRFaces(X1+dx,X2,tol);
    mesh.tagFacesInDomain(fRFaces,fDeadZone,RIGHT_DEAD_ZONE);

  }
  
}


 void LaxFriedrichsSystemMPI::iterate(double t,double tEnd) 
{
  ConservativeMethodForSystem::iterate(t,tEnd);
  updateDeadZone();
}

void LaxFriedrichsSystemMPI::updateGhostCells() 
{
  for (unsigned i=0;i<n_components;i++)
  {
    NetMPI::exchangeDataRedBlack(LGhostCellSndMap,LGhostCellSnd[i],LGhostCellRcvMap,LGhostCellRcv[i],
                                 RGhostCellSndMap,RGhostCellSnd[i],RGhostCellRcvMap,RGhostCellRcv[i],
                                 m_sol[i],SEND_GHOSTCELL);
  }
}




void LaxFriedrichsSystemMPI::updateDeadZone() 
{
  for (unsigned i=0;i<n_components;i++)
  {
    NetMPI::exchangeDataRedBlack(cLDeadZoneSndMap,cLDeadZoneSnd[i],cLDeadZoneRcvMap,cLDeadZoneRcv[i],
                                 cRDeadZoneSndMap,cRDeadZoneSnd[i],cRDeadZoneRcvMap,cRDeadZoneRcv[i],
                                 m_sol[i],DEAD_ZONE_VALUE);
  }
}


void LaxFriedrichsSystemMPI::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  vQPos(n_components),vQNeg(n_components),vBC(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();
    unsigned nFaces = m_mesh.numFaces();

    //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
    VecTag::iterator itDZ = fDeadZone.begin();
    for (unsigned faceIndex=0;faceIndex < nFaces; faceIndex++,face++,itDZ++)
    {
      if (*itDZ != ACTIVE_ZONE)
        continue;

      //Get the values of the previous solution in the right and left of the face
      this->getValuesOfTheFaceCells(m_mesh,m_fBC,m_solAnt,face,vQNeg,vQPos);

      //Get the boundary conditions in that face.
      //If the face doesnt have boundary condition for the ith
      //componente, then vBC == infinity.
      getBC(faceIndex,vBC);

      
      //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
      face->getValidCellIndices(posCell,negCell);

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

      //Now get the flux component by component
      for (unsigned cmp=0;cmp<n_components;cmp++)
      {
        
        //See if the face has prescribed boundary condition BC      
        double flux;
        if (vBC(cmp) == INFINITY)    
        {
          /*
            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() * m_flux.fluxAtFace(face,faceIndex,posCell,negCell,vQPos,vQNeg,cmp) - mPor*(vQNeg(cmp)-vQPos(cmp))/6.0;
        }
        else //the face has BC, so dont add the stability term
        {
          flux = dt*face->areaPerCellVol() * m_flux.fluxAtFace(face,faceIndex,posCell,negCell,vBC,vBC,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[cmp](face->getPosCell()) +=  - flux/posPor;

        if (face->hasNegCell())
          m_sol[cmp](face->getNegCell()) +=  + flux/negPor;
        
      } //end for each component
    } //end for each face
    updateGhostCells();
  } //end for each time step
}


 void LaxFriedrichsSystemMPI::updateDataForDynamicModule() 
{
  updateDeadZone();
}


 double LaxFriedrichsSystemMPI::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 NetMPI::getMinValue(dt);
}

---