laxfriedrichsmethodmpi.cpp

Go to the documentation of this file.

#include "laxfriedrichsmethodmpi.h"
#include "netmpi.h"
#include "fxyslabid.h"
#include "numericmethods.h"

#define SEND_GHOSTCELL 201
#define DEAD_ZONE_VALUE 202


void LaxFriedrichsMethodMPI::setCommunicationBuffersWithNoOverlap(OrthoMesh &mesh) 
{
  unsigned elemsInXY = mesh.numElemY()*mesh.numElemZ();
  double dx = mesh.get_dx();
  double tol = dx/4.0;

  //If we dont have mesh overlap. each subdomain gets boundary conditions from your neighbors
  //and there is no dead zones
  if (NetMPI::rank() != 0)
  {
    double X0 = mesh.getP()[0];


    //With no dead zones each subdomain send the cells at the shared boundary
    LGhostCellSndMap.reserve(elemsInXY);
    FXYSlabId fSnd(X0,X0+dx,tol);
    mesh.getCellsIndicesInDomain(fSnd,LGhostCellSndMap);

    //Now each subdomain receive cells values coming from the other side of its boundaries.
    //These values were sent by its neighbours domains.
    LGhostCellRcvMap.reserve(elemsInXY);
    FXYSlabId fRcv(X0,X0,tol);
    mesh.getFacesIndicesInDomain(fRcv,LGhostCellRcvMap);

    //Initialize buffers to send and get data.
    LGhostCellRcv.reinit(LGhostCellRcvMap.size());
    LGhostCellSnd.reinit(LGhostCellSndMap.size());
  }
  if (!NetMPI::isLastProcess())
  {
    double X0 = mesh.getQ()[0];

    //With no dead zones each subdomain send the cells at the shared boundary
    RGhostCellSndMap.reserve(elemsInXY);
    FXYSlabId fSnd(X0,X0-dx,tol);
    mesh.getCellsIndicesInDomain(fSnd,RGhostCellSndMap);

    //Now each subdomain receive cells values coming from the other side of its boundaries.
    //These values were sent by its neighbours domains.
    RGhostCellRcvMap.reserve(elemsInXY);
    FXYSlabId fRcv(X0,X0,tol);
    mesh.getFacesIndicesInDomain(fRcv,RGhostCellRcvMap);

    //Initialize buffers to send and get data.
    RGhostCellRcv.reinit(RGhostCellRcvMap.size());
    RGhostCellSnd.reinit(RGhostCellSndMap.size());
  }
}



void LaxFriedrichsMethodMPI::setCommunicationBuffers(OrthoMesh &mesh,int mesh_overlap) 
{
  if (mesh_overlap == 0)
    return setCommunicationBuffersWithNoOverlap(mesh);

  //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);
      cLDeadZoneRcv.reinit(cLDeadZoneRcvMap.size());
      cLDeadZoneSnd.reinit(cLDeadZoneSndMap.size());

      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);
      LGhostCellRcv.reinit(LGhostCellRcvMap.size());
      LGhostCellSnd.reinit(LGhostCellSndMap.size());
      assert(LGhostCellRcvMap.size() == elemsInXY);
      assert(LGhostCellSndMap.size() == elemsInXY);

      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);
      cRDeadZoneRcv.reinit(cRDeadZoneRcvMap.size());
      cRDeadZoneSnd.reinit(cRDeadZoneSndMap.size());

      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);
      RGhostCellRcv.reinit(RGhostCellRcvMap.size());
      RGhostCellSnd.reinit(RGhostCellSndMap.size());
      assert(RGhostCellRcvMap.size() == elemsInXY);
      assert(RGhostCellSndMap.size() == elemsInXY);

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

  











 LaxFriedrichsMethodMPI::~LaxFriedrichsMethodMPI() 
{

}

 LaxFriedrichsMethodMPI::LaxFriedrichsMethodMPI(OrthoMesh &mesh,Function3D &fInitU,const VecDouble &cPor,Function3D &fPrescribedU,FaceFluxFunction &flux, double CFL,int mesh_overlap)
   :LaxFriedrichsMethod(mesh,fInitU,cPor,fPrescribedU,flux,CFL),mesh_overlap(mesh_overlap)
 {
   NetMPI::trace("Entering Communication Buffers ");
   setCommunicationBuffers(mesh,mesh_overlap);
   NetMPI::trace("End");
   
}


 void LaxFriedrichsMethodMPI::iterate(double t, double tEnd) 
{
  unsigned posCell,negCell;
  VecReq _req(4);
  double timeInterval = tEnd-t;
  double localDt = this->calculateTimeStep(m_mesh,timeInterval,m_CFL,getPorosity(),m_flux);
  double dt=negotiateTimeStep(localDt);
  NetMPI::trace("Lax Negotiating dt: Sent %g, got %g",localDt,dt); 
  
  //Get the number of iterations
  int nIteraction = (int) round(timeInterval/dt); //number of iteractions

  //For each time step
  for (int count = 0; count < nIteraction; count++)
  {
    //Get the face iterator and the number of faces.
    OrthoMesh::Face_It face = m_mesh.begin_face();
    unsigned nFaces = m_mesh.numFaces();
    
    m_cValuesPrev=m_cValues;

    //For each face
    VecTag::iterator itDZ = fDeadZone.begin();
    for (unsigned faceIndex=0;faceIndex < nFaces; faceIndex++,face++,itDZ++)
    {
      //Check if the face is in active zone
      if (*itDZ != ACTIVE_ZONE)
        continue;
      
      //Get the values of the previous solution in the right and left of the face
      double SwPos,SwNeg;
      this->getValuesOfTheFaceCells(m_mesh,m_fValues,m_cValuesPrev,face,SwNeg,SwPos);

      //Get the values 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;
      
      //See if the face has prescribed boundary condition BC      
      double flux;
      double bc = (getBC(faceIndex));
      if (bc == 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,SwPos,SwNeg) - mPor*(SwNeg-SwPos)/6.0;
      }
      else //the face has BC, so dont add the stability term
      {
        flux = dt*face->areaPerCellVol() * m_flux.fluxAtFace(face,faceIndex,posCell,negCell,bc,bc);
      }

      
      
      //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_cValues(face->getPosCell()) +=  - flux/posPor;

      if (face->hasNegCell())
        m_cValues(face->getNegCell()) +=  + flux/negPor;
    }
    updateGhostCells();
    
  }
  updateDeadZone();
}







void LaxFriedrichsMethodMPI::updateGhostCells() 
{
  if (mesh_overlap)
    NetMPI::exchangeDataRedBlack(LGhostCellSndMap,LGhostCellSnd,LGhostCellRcvMap,LGhostCellRcv,
                                 RGhostCellSndMap,RGhostCellSnd,RGhostCellRcvMap,RGhostCellRcv,
                                 m_cValues,SEND_GHOSTCELL);
  else
    NetMPI::exchangeDataRedBlack(LGhostCellSndMap,LGhostCellSnd,LGhostCellRcvMap,LGhostCellRcv,
                                 RGhostCellSndMap,RGhostCellSnd,RGhostCellRcvMap,RGhostCellRcv,
                                 m_cValues,m_fValues,SEND_GHOSTCELL);
    


}

void LaxFriedrichsMethodMPI::updateDeadZone() 
{
  if (mesh_overlap)//If no mesh_overlap, we dont have dead zones
  {
    NetMPI::exchangeDataRedBlack(cLDeadZoneSndMap,cLDeadZoneSnd,cLDeadZoneRcvMap,cLDeadZoneRcv,
                                 cRDeadZoneSndMap,cRDeadZoneSnd,cRDeadZoneRcvMap,cRDeadZoneRcv,
                                 m_cValues,DEAD_ZONE_VALUE);
  }
}


void LaxFriedrichsMethodMPI::unitTest() 
{
  std::ostream &out = NetMPI::getLog();
  out << "\n\nLeft Ghost Rcv\n";
  m_mesh.printCells(LGhostCellRcvMap,out);
  out << "\n\nLeft Ghost Snd\n";
  m_mesh.printCells(LGhostCellSndMap,out);
  out << "\n\nRight Ghost Rcv\n";
  m_mesh.printCells(RGhostCellRcvMap,out);
  out << "\n\nRight Ghost Snd\n";
  m_mesh.printCells(RGhostCellSndMap,out);

}


double LaxFriedrichsMethodMPI::negotiateTimeStep(double dt) 
{
  double dd;
  MPI_Allreduce(&dt,&dd,1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD);
  return dd;
}

---