hybridpressionmodulempi.cpp

Go to the documentation of this file.

#include "hybridpressionmodulempi.h"
#include "netmpi.h"
#include "transportbase.h"
#include "numericmethods.h"
#include "hdf5dealwriter.h"
#include "hdf5orthowriter.h"

#define BLACK_TO_RED 103
#define RED_TO_BLACK 104
#define B_PARAMETER 101
#define BC_MESSAGE 102
#define ROBIN_BC_MESSAGE 107


HybridPressionModuleMPI::HybridPressionModuleMPI(OrthoMesh &mesh,Function3D &fPrescribedVelocities, Function3D &fPrescribedPression, Function3D &K,Function1D &fMobT,FunctionOfCellFields &fGravSource,double Xsize,double b,double theta,int meshOverlap,double maxTol, int debugLevel,std::ostream &out) 
  :HybridPressionModule(mesh,fPrescribedVelocities,fPrescribedPression,K,fMobT,fGravSource,debugLevel,false),
   m_out(out)
{
  assert(m_theta >= 0.0 && m_theta <=1.0);
  m_max_tol = maxTol;
  m_meshOverlap = meshOverlap;
  setCommunication();
  m_xSize = Xsize;
  m_b=b;
  m_sol=0.0;
  m_solPrev=m_sol;
  m_theta=theta;
  if (meshOverlap<=0)
    throw new Exception("Invalid meshOverlap parameter");

#ifdef DEBUG
  checkBC(fPrescribedVelocities,fPrescribedPression);
#endif

  //Retrieve this code in the final release
    VecDouble Vv(m_mesh.numVertices());
    HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 

    
}

void HybridPressionModuleMPI::setCommunication() 
{
  if (NetMPI::rank()!=0)
  {
    m_mesh.getFacesInYZSlab(m_mesh.getP()[0],fLeftBCRcvMap);
    m_mesh.getFacesInYZSlab(m_mesh.getP()[0]+2*m_meshOverlap*m_mesh.getDX()[0],fLeftBCSndMap);

    //get B parameter to send to the left domain
    BL_Snd.reinit(fLeftBCRcvMap.size());
    BL_Rcv.reinit(fLeftBCRcvMap.size());
    LBC_Snd.reinit(fLeftBCRcvMap.size());
    LBC_Rcv.reinit(fLeftBCRcvMap.size());
    LBC_Prev.reinit(fLeftBCRcvMap.size());

    LBC_Rcv=0.0;
    LBC_Prev=0.0;
    
  }
  if (!NetMPI::isLastProcess())
  {
    m_mesh.getFacesInYZSlab(m_mesh.getQ()[0],fRightBCRcvMap);
    m_mesh.getFacesInYZSlab(m_mesh.getQ()[0]-2*m_meshOverlap*m_mesh.getDX()[0],fRightBCSndMap);

    //get B parameter to send to the right domain
    BR_Snd.reinit(fRightBCRcvMap.size());
    BR_Rcv.reinit(fRightBCRcvMap.size());
    RBC_Snd.reinit(fRightBCRcvMap.size());
    RBC_Rcv.reinit(fRightBCRcvMap.size());
    RBC_Prev.reinit(fRightBCRcvMap.size());
    RBC_Rcv=0.0;
    RBC_Prev=0.0;
  }
  if (NetMPI::rank() == 0)
  {
    m_vTol.reinit(NetMPI::nProcess());
  }
  
  

}




void HybridPressionModuleMPI::checkBC(Function3D &fPrescribedVelocities, Function3D &fPrescribedPression) 
{
  Point3D p;  

  OrthoMesh::Cell_It cell = m_mesh.begin_cell(),
    endc = m_mesh.end_cell();

  
   //For eahc cell in the Mesh
   for(;cell != endc;cell++)
   {
     //For each boundary face of the cell
     for (unsigned int face_no = 0; face_no < GeometryInfo<DIM>::faces_per_cell;++face_no)
     {
       FaceDirection3D fDir = (FaceDirection3D) face_no;
       OrthoMesh::Face_It face = cell->face(fDir);
       if (face->at_boundary())
       {
         bool found=false;
         face->barycenter(p);
         if (fPrescribedPression.isInDomain(p))
         {
           found=true;
         }
         if (fPrescribedVelocities.isInDomain(p))
         {
           if (found)
             throw new Exception("Boundary conditions clash");
           else
             found=true;
         }
         if ( (std::find(fLeftBCRcvMap.begin(),fLeftBCRcvMap.end(),face->index()) !=fLeftBCRcvMap.end()) ||
              (std::find(fRightBCRcvMap.begin(),fRightBCRcvMap.end(),face->index()) !=fRightBCRcvMap.end()))
         {
           if (found)
             throw new Exception("Boundary conditions clash");
           else
             found=true;
         }
         if (!found)
         {
           Point3D p;
           face->barycenter(p); 
           throw new Exception("Found face I %d N %u (%u,%u,%u)-%g,%g,%g without boundary condition",face->index(),static_cast<unsigned>(face->getNormalOrientation()),face->getI(),face->getJ(),face->getK(),p[0],p[1],p[2]);
         }
           
       }
     }
   }
}


void HybridPressionModuleMPI::iterateRed(TransportBase &trans) 
{
  //SET AND GET  MPI VARIABLES 
  MPI_Request _req[4];
  MPI_Status stats[4];
  _req[0]=_req[1]=_req[2]=_req[3]=MPI_REQUEST_NULL;

  double tol;
  int rank = NetMPI::rank();
  int next = rank+1;
  int prev = rank-1;


  const VecDouble &S1 = trans.getSolutionAtCells(0);

  assemblySystem(trans,m_A,m_B,m_K,S1,m_fMobT,m_fMobGrav);
  initIteration();
  

  m_itCnt=0;
  do {
    m_RHS = m_B;
    addRobinBC(LBC_Rcv,BL_Rcv,fLeftBCRcvMap,m_K,S1,m_fMobT,m_RHS);
    addRobinBC(RBC_Rcv,BR_Rcv,fRightBCRcvMap,m_K,S1,m_fMobT,m_RHS);
    m_solPrev=m_sol;
    solve();

    //Send left and right robin condition
    getRobinBC(LBC_Snd,m_K,m_sol,S1,m_fMobT,fLeftBCSndMap,true);
    NetMPI::Isend(LBC_Snd,prev,ROBIN_BC_MESSAGE,_req[0]);
    getRobinBC(RBC_Snd,m_K,m_sol,S1,m_fMobT,fRightBCSndMap,false);
    NetMPI::Isend(RBC_Snd,next,ROBIN_BC_MESSAGE,_req[1]);

    //Receive left and right  robin condition
    RBC_Prev.swap(RBC_Rcv);
    LBC_Prev.swap(LBC_Rcv);
    NetMPI::Irecv(LBC_Rcv,prev,ROBIN_BC_MESSAGE,_req[2]);
    NetMPI::Irecv(RBC_Rcv,next,ROBIN_BC_MESSAGE,_req[3]);
    MPI_Waitall(4,_req,stats);
 
    if (rank != 0)
      LBC_Rcv.sadd((1.0-m_theta),m_theta,LBC_Prev);
    if (!NetMPI::isLastProcess())
      RBC_Rcv.sadd((1.0-m_theta),m_theta,RBC_Prev);

    //    tol = std::max(NumericMethods::relError(LBC_Rcv,LBC_Prev),
    //             NumericMethods::relError(RBC_Rcv,RBC_Prev));
    tol = NumericMethods::relError(m_solPrev,m_sol);
  } while (checkConvergence(tol));

  getNormalVelocitiesFromPressure(m_Vq,m_K,m_sol,S1,m_fMobT,m_fMobGrav,LBC_Prev,RBC_Prev);

}



void HybridPressionModuleMPI::iterateBlack(TransportBase &trans) 
{
  //SET AND GET  MPI VARIABLES 
  MPI_Request _req[4];
  MPI_Status stats[4];
  _req[0]=_req[1]=_req[2]=_req[3]=MPI_REQUEST_NULL;

  double tol;
  int rank = NetMPI::rank();
  int next = rank+1;
  int prev = rank-1;

  const VecDouble &S1 = trans.getSolutionAtCells(0);
  assemblySystem(trans,m_A,m_B,m_K,S1,m_fMobT,m_fMobGrav);
  initIteration();
  do {

    m_RHS = m_B;
    addRobinBC(LBC_Rcv,BL_Rcv,fLeftBCRcvMap,m_K,S1,m_fMobT,m_RHS);
    addRobinBC(RBC_Rcv,BR_Rcv,fRightBCRcvMap,m_K,S1,m_fMobT,m_RHS);
    m_solPrev=m_sol;
    solve();

    //Receive right and left  robin condition
    RBC_Prev.swap(RBC_Rcv);
    LBC_Prev.swap(LBC_Rcv);
    NetMPI::Irecv(RBC_Rcv,next,ROBIN_BC_MESSAGE,_req[0]);
    NetMPI::Irecv(LBC_Rcv,prev,ROBIN_BC_MESSAGE,_req[1]);

    //Send right and left robin condition
    getRobinBC(RBC_Snd,m_K,m_sol,S1,m_fMobT,fRightBCSndMap,false);
    NetMPI::Isend(RBC_Snd,next,ROBIN_BC_MESSAGE,_req[3]);
    getRobinBC(LBC_Snd,m_K,m_sol,S1,m_fMobT,fLeftBCSndMap,true);
    NetMPI::Isend(LBC_Snd,prev,ROBIN_BC_MESSAGE,_req[4]);
    MPI_Waitall(4,_req,stats);

    
    if (rank != 0)
      LBC_Rcv.sadd((1.0-m_theta),m_theta,LBC_Prev);
    if (!NetMPI::isLastProcess())
      RBC_Rcv.sadd((1.0-m_theta),m_theta,RBC_Prev);

    //tol = std::max(NumericMethods::relError(LBC_Rcv,LBC_Prev),
    //             NumericMethods::relError(RBC_Rcv,RBC_Prev));

    tol = NumericMethods::relError(m_solPrev,m_sol);


  }while (checkConvergence(tol));

  getNormalVelocitiesFromPressure(m_Vq,m_K,m_sol,S1,m_fMobT,m_fMobGrav,LBC_Rcv,RBC_Rcv);
  
}



 


  
  //is it a black node ?
  



 HybridPressionModuleMPI::~HybridPressionModuleMPI() 
{

}

void HybridPressionModuleMPI::exchangeData(VecDouble &lSnd,VecDouble &lRcv,VecDouble &rSnd,VecDouble &rRcv) 
{


  int rank = NetMPI::rank();
  int next=rank+1;
  int prev=rank-1;
  MPI_Request reqs[4];
  reqs[0]=reqs[1]=reqs[2]=reqs[3]=MPI_REQUEST_NULL;
  
  MPI_Status stats;

  assert(!(rank == 0 && lSnd.size() != 0));
  assert(!(NetMPI::isLastProcess() && rSnd.size() != 0));
  
  if (rank%2 == 0) //black domain
  {
    if (rank!= 0)
      MPI_Send(&(lSnd(0)),lSnd.size(),MPI_DOUBLE,prev,BLACK_TO_RED,MPI_COMM_WORLD);
    if (!NetMPI::isLastProcess())
      MPI_Send(&(rSnd(0)),rSnd.size(),MPI_DOUBLE,next,BLACK_TO_RED,MPI_COMM_WORLD);

    if (rank!=0)
      MPI_Recv(&(lRcv(0)),lRcv.size(),MPI_DOUBLE,prev,RED_TO_BLACK,MPI_COMM_WORLD,&stats);
    if (!NetMPI::isLastProcess())
      MPI_Recv(&(rRcv(0)),rRcv.size(),MPI_DOUBLE,next,RED_TO_BLACK,MPI_COMM_WORLD,&stats);
  }
  else
  {
    if (!NetMPI::isLastProcess())
      MPI_Recv(&(rRcv(0)),rRcv.size(),MPI_DOUBLE,next,BLACK_TO_RED,MPI_COMM_WORLD,&stats);
    MPI_Recv(&(lRcv(0)),lRcv.size(),MPI_DOUBLE,prev,BLACK_TO_RED,MPI_COMM_WORLD,&stats);

    if (!NetMPI::isLastProcess())
      MPI_Send(&(rSnd(0)),rSnd.size(),MPI_DOUBLE,next,BLACK_TO_RED,MPI_COMM_WORLD);
    MPI_Send(&(lSnd(0)),lSnd.size(),MPI_DOUBLE,prev,BLACK_TO_RED,MPI_COMM_WORLD);

    

  }

}


void HybridPressionModuleMPI::getBParameter(VecDouble &vB, const VecDouble &K,const VecDouble &S1,Function1D &fMobT,const VecIndex &fI) 
{
  assert(vB.size() == fI.size());
  for (unsigned i=0;i<vB.size();i++)
  {
    OrthoMesh::Face_It face = m_mesh.get_face(fI[i]);
    assert(!face->at_boundary());
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);

    double K1 = K(c1)*fMobT(S1(c1));
    double K2 = K(c2)*fMobT(S1(c2));
     
    double Keff = 2*K1*K2/(K1 + K2);
    vB(i)=m_b*m_xSize/Keff;
  }

}


void HybridPressionModuleMPI::getRobinBC(VecDouble &vRobin,const VecDouble &K,const VecDouble &vP,const VecDouble &S1,Function1D &fMobT,const VecIndex &fI, bool bNegateB) 
{
  assert(vRobin.size() == fI.size());
  double b = bNegateB ? -m_b*m_xSize : m_b*m_xSize;
  for (unsigned i=0;i<vRobin.size();i++)
  {
    OrthoMesh::Face_It face = m_mesh.get_face(fI[i]);
    assert(!face->at_boundary());
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);
    double p1 = vP(c1);
    double p2 = vP(c2);
    double K1 = K(c1)*fMobT(S1(c1));
    double K2 = K(c2)*fMobT(S1(c2));
    vRobin(i) = +b*(p2-p1)*face->areaPerCellVol() + (K1*p1 + K2*p2)/(K1 + K2);
  }
}


 void HybridPressionModuleMPI::addRobinBC(const VecDouble &v1,const VecDouble &vB,const VecIndex &fI, const VecDouble &vK,const VecDouble &S1,Function1D &fMobT,VecDouble &vRHS) 
{
  assert(vB.size() == fI.size());
  assert(v1.size() == fI.size());
  assert(vK.size() == m_mesh.numCells());
  assert(S1.size() == m_mesh.numCells());
  assert(vRHS.size() == m_mesh.numCells());
 
  for (unsigned i=0;i<vB.size();i++)
  {
    OrthoMesh::Face_It face = m_mesh.get_face(fI[i]);
    if (!face->at_boundary())
    {
      std::cout << "throw exception\n" << std::endl;
      throw new Exception("Problem at face %d",fI[i]);

    }
    assert(face->at_boundary());
    assert(face->getNormalOrientation() == OrthoMesh::NORMAL_X);
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);
    double fArea = face->area(); 
    if (c2 != OrthoMesh::INVALID_INDEX)
    {
      double K = vK(c2)*fMobT(S1(c2));
      vRHS(c2)+= K*2*v1(i)*fArea/(m_mesh.get_dx() + 2*K*vB(i));
    }
    if (c1 != OrthoMesh::INVALID_INDEX)
    {
      double K = vK(c1)*fMobT(S1(c1));
      vRHS(c1)+= K*2*v1(i)*fArea/(m_mesh.get_dx() + 2*K*vB(i));
    }
    
    
  }
}


void HybridPressionModuleMPI::addRobinBC(const VecDouble &vB,const VecIndex &fI, const VecDouble &vK,const VecDouble &S1,Function1D &fMobT,SparseMatrix<double> &A) 
{
  assert(vB.size() == fI.size());
  assert(vK.size() == m_mesh.numCells());
  assert(S1.size() == m_mesh.numCells());
  assert(A.m() == A.n() && A.m() == m_mesh.numCells());

  
  for (unsigned i=0;i<vB.size();i++)
  {
    OrthoMesh::Face_It face = m_mesh.get_face(fI[i]);
    assert(face->at_boundary());
    assert(face->getNormalOrientation() == OrthoMesh::NORMAL_X);
    unsigned c1,c2;
    face->getAdjCellIndices(c1,c2);
    double fArea = face->area(); 
    if (c2 != OrthoMesh::INVALID_INDEX) //left boundary
    {
      double K = vK(c2)*fMobT(S1(c2));
      double BC =  2*K*fArea/(m_mesh.get_dx() + 2*K*vB(i));
      A.add(c2,c2,BC);
    }
    if (c1 != OrthoMesh::INVALID_INDEX)//right boundary
    {
      double K = vK(c1)*fMobT(S1(c1));
      double BC = 2*K*fArea/(m_mesh.get_dx() + 2*K*vB(i));
      A.add(c1,c1,BC);
    }
    
    
  }
  
}


 void HybridPressionModuleMPI::iterate(TransportBase &trans) 
{
  int rank = NetMPI::rank();

  if (rank%2 == 0)
  {
    iterateRed(trans);

  }
  else if (rank%2 == 1)
  {
    iterateBlack(trans);
  }
    
}




void HybridPressionModuleMPI::printOutput() 
{
  if (_debugLevel > 0)
  {
    VecDouble v(m_mesh.numCells());
    VecDouble Vv(m_mesh.numVertices());
    HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 

    this->getPressionAtCells(v);
    m_mesh.projectCentralValuesAtVertices(v,Vv);
    hdf5.writeScalarField(Vv,"P");

    //    GnuPlotAnim &plot = GnuPlotAnim::getGnuPlotAnim();
    //plot.plotVerticeValuesAtFixedZ(m_mesh.getTriangulation(),Vv,"pression", "with lines",250);
    //plot.NextScene();

  
    if (_debugLevel > 1)
    {
      Matrix Mvel(m_mesh.numCells(),3);
      this->getVelocitiesAtCells(Mvel);
      if (_debugLevel > 2)
      {
        printf("\nVelocities at cells' barycenters\n\n");
        Mvel.print_formatted(m_out,3,true,0,"0");
      }
      VecDouble v(m_mesh.numCells());
      for (unsigned i=0;i<v.size();i++)
      {
        v(i) = Mvel(i,0); 
      }
      VecDouble minVel(3);
      VecDouble maxVel(3);
      NumericMethods::getMinMaxValueByColumn(Mvel,minVel,maxVel);
      printf("Min Vel <%g,%g,%g> - Max Vel <%g,%g,%g>\n",minVel(0),minVel(1),minVel(2),maxVel(0),maxVel(1),maxVel(2));
  
      m_mesh.projectCentralValuesAtVertices(v,Vv);
      hdf5.writeScalarField(Vv,"Vel");
    }
  }

}


void HybridPressionModuleMPI::getNormalVelocitiesFromPressure(VecDouble &Vq,const VecDouble &K,const VecDouble &vP,const VecDouble &S1,Function1D &fMobT,FunctionOfCellFields &fMobGrav,const VecDouble &LRobinBC,const VecDouble &RRobinBC) 
{
  HybridPressionModule::getNormalVelocitiesFromPressure(Vq,K,vP,S1,fMobT,fMobGrav);
  assert(RRobinBC.size() ==0 || RRobinBC.size() == fRightBCRcvMap.size());
  assert(LRobinBC.size() ==0 || LRobinBC.size() == fLeftBCRcvMap.size());

  for (unsigned i=0;i<LRobinBC.size();i++)
  {
    int faceIndex = fLeftBCRcvMap[i];
    OrthoMesh::Face_It face = m_mesh.get_face(faceIndex);
    unsigned  c2 = face->getNegCell();
    assert(face->at_boundary());
    assert(face->getNormalOrientation() == OrthoMesh::NORMAL_X);
    assert(c2!=OrthoMesh::INVALID_INDEX);
    double k = K(c2)*fMobT(S1(c2));
    double p = vP(c2);
    Vq(faceIndex) = - 2*k*(p-LRobinBC(i))/(m_mesh.get_dx() + 2*k*BL_Rcv(i));
  }
  
  for (unsigned i=0;i<RRobinBC.size();i++)
  {
    int faceIndex = fRightBCRcvMap[i];
    OrthoMesh::Face_It face = m_mesh.get_face(faceIndex);
    unsigned c1 = face->getPosCell();
    assert(face->at_boundary());
    assert(face->getNormalOrientation() == OrthoMesh::NORMAL_X);
    assert(c1!=OrthoMesh::INVALID_INDEX);
    double k = K(c1)*fMobT(S1(c1));
    double p = vP(c1);
    Vq(faceIndex) = - 2*k*(RRobinBC(i)-p)/(m_mesh.get_dx() + 2*k*BR_Rcv(i));
  }
}


void HybridPressionModuleMPI::assemblySystem(TransportBase &trans,SparseMatrix<double> &A,VecDouble &B,VecDouble &K,const VecDouble &S1, Function1D &fMobT, FunctionOfCellFields &fMobGrav) 
{
  //SET AND GET  MPI VARIABLES 
  MPI_Request req[4];
  MPI_Status stats[4];
  req[0]=req[1]=req[2]=req[3]=MPI_REQUEST_NULL;
   
  int rank = NetMPI::rank();
  int next = rank+1;
  int prev = rank-1;

  
  if (rank%2 == 0)
  {

    //First of all we need to get the B parameter used in the domain decomposition

    //Send B to left and right domain
    getBParameter(BL_Snd,K,S1,fMobT,fLeftBCSndMap);
    NetMPI::Isend(BL_Snd,prev,B_PARAMETER,req[0]);
    getBParameter(BR_Snd,K,S1,fMobT,fRightBCSndMap);
    NetMPI::Isend(BR_Snd,next,B_PARAMETER,req[1]);

    //Receive B from left  and right  domain
    NetMPI::Irecv(BL_Rcv,prev,B_PARAMETER,req[2]);
    NetMPI::Irecv(BR_Rcv,next,B_PARAMETER,req[3]);

    //Set the gravity mobility function defined in the cells
    //with the solutions of the transport module
    for (unsigned fieldId=0;fieldId<fMobGrav.numFields();fieldId++) //Set the fields of the gravity source term
      fMobGrav.setField(trans.getSolutionAtCells(fieldId),fieldId);

    //assembly system
    HybridPressionModule::assemblySystem(A,B,K,S1,fMobT,fMobGrav);


    //Sync point reached
    //after this point all processors have the B parameter from its neighbors
    MPI_Waitall(4,req,stats);

    addRobinBC(BL_Rcv,fLeftBCRcvMap,K,S1,fMobT,A);
    addRobinBC(BR_Rcv,fRightBCRcvMap,K,S1,fMobT,A);

  
  }
  else if (rank%2 == 1)
  {

    //First of all we need to get the B parameter used in the domain decomposition

    //Receive B from right and left domain
    NetMPI::Irecv(BR_Rcv,next,B_PARAMETER,req[0]);
    NetMPI::Irecv(BL_Rcv,prev,B_PARAMETER,req[1]);

    //Send B to right and left domain
    getBParameter(BR_Snd,K,S1,fMobT,fRightBCSndMap);
    NetMPI::Isend(BR_Snd,next,B_PARAMETER,req[3]);
    getBParameter(BL_Snd,K,S1,fMobT,fLeftBCSndMap);
    NetMPI::Isend(BL_Snd,prev,B_PARAMETER,req[2]);

    //Set the gravity mobility function defined in the cells
    //with the solutions of the transport module
    for (unsigned fieldId=0;fieldId<fMobGrav.numFields();fieldId++) //Set the fields of the gravity source term
      fMobGrav.setField(trans.getSolutionAtCells(fieldId),fieldId);

    //assembly system
    HybridPressionModule::assemblySystem(A,B,K,S1,fMobT,fMobGrav);


    //Sync point reached
    //after this point all processors have the B parameter from its neighbors
    MPI_Waitall(4,req,stats);

    //With the B parameter we can add the robin conditions to the system matrix.
    addRobinBC(BL_Rcv,fLeftBCRcvMap,K,S1,fMobT,A);
    addRobinBC(BR_Rcv,fRightBCRcvMap,K,S1,fMobT,A);
  }
}

int HybridPressionModuleMPI::checkConvergence(double tol) 
{
  int bContinue;
  HDF5OrthoWriter::getHDF5OrthoWriter().setVariable("ItCnt",static_cast<double>(m_itCnt));

  if (NetMPI::rank() != 0)
  {
    MPI_Gather(&tol,1,MPI_DOUBLE,NULL,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
    MPI_Bcast(&bContinue,1,MPI_INT,0,MPI_COMM_WORLD);
  }
  else
  {
    m_vTol(0)=tol;
    MPI_Gather(&tol,1,MPI_DOUBLE,&(m_vTol(0)),1,MPI_DOUBLE,0,MPI_COMM_WORLD);
    bContinue=false;
    printf("%d Tolerance: ",m_itCnt++);
    for (unsigned i=0;i<m_vTol.size();i++)
    {
      printf("%g ",m_vTol(i));
      if (m_vTol(i) >= m_max_tol)
      {
        bContinue = true;
        break;
      }
    }
    printf("\n");
    MPI_Bcast(&bContinue,1,MPI_INT,0,MPI_COMM_WORLD);
  }
  return bContinue;
}

---