HybridPressionModuleMPI Class Reference

#include <hybridpressionmodulempi.h>

Inheritance diagram for HybridPressionModuleMPI:
Inheritance graph
[legend]
Collaboration diagram for HybridPressionModuleMPI:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 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=std::cout)
virtual ~HybridPressionModuleMPI ()
virtual void iterate (TransportBase &trans)
virtual void printOutput ()

Protected Member Functions

void setCommunication ()
void checkBC (Function3D &fPrescribedVelocities, Function3D &fPrescribedPression)
void exchangeData (VecDouble &lSnd, VecDouble &lRcv, VecDouble &rSnd, VecDouble &rRcv)
void getBParameter (VecDouble &vB, const VecDouble &K, const VecDouble &S1, Function1D &fMobT, const VecIndex &fI)
void addRobinBC (const VecDouble &v1, const VecDouble &pB, const VecIndex &fI, const VecDouble &vK, const VecDouble &S1, Function1D &fMobT, VecDouble &vRHS)
void addRobinBC (const VecDouble &pB, const VecIndex &fI, const VecDouble &vK, const VecDouble &S1, Function1D &fMobT, SparseMatrix< double > &A)
void getRobinBC (VecDouble &vRobin, const VecDouble &K, const VecDouble &vP, const VecDouble &S1, Function1D &fMobT, const VecIndex &fI, bool bNegateB)
virtual void iterateRed (TransportBase &trans)
virtual void iterateBlack (TransportBase &trans)
void assemblySystem (TransportBase &trans, SparseMatrix< double > &A, VecDouble &B, VecDouble &K, const VecDouble &S1, Function1D &fMobT, FunctionOfCellFields &fMobGrav)
int checkConvergence (double tol)
void getNormalVelocitiesFromPressure (VecDouble &Vq, const VecDouble &K, const VecDouble &vP, const VecDouble &S1, Function1D &fMobT, FunctionOfCellFields &fMobGrav, const VecDouble &LRobinBC, const VecDouble &RRobinBC)
virtual void initIteration ()=0
virtual void solve ()=0

Protected Attributes

std::ostream & m_out
VecIndex fLeftBCRcvMap
VecIndex fRightBCRcvMap
VecIndex fLeftBCSndMap
VecIndex fRightBCSndMap
int m_meshOverlap
double m_xSize
double m_b
double m_theta
unsigned m_itCnt
VecDouble m_RHS
VecDouble m_solPrev
VecDouble BL_Snd
VecDouble BL_Rcv
VecDouble BR_Snd
VecDouble BR_Rcv
VecDouble RBC_Snd
VecDouble RBC_Rcv
VecDouble LBC_Snd
VecDouble LBC_Rcv
VecDouble RBC_Prev
VecDouble LBC_Prev
VecDouble m_vTol
double m_max_tol

Detailed Description

This class implements the hybrid solver with domain decomposition. The domain is splitted with adjacent subdomains along X direction ___________________ |R|B|R|B.... | -------------------

Adjacent domains have a non zero intersection containing at least one Slab of elements in the YZ plane. Processes change robin boundary conditions along the interfaces and compute their local problems until they reach a convergence. More information consult the Doctor's thesis of Oscar Osmay Delgado, UERJ, 2001

Definition at line 18 of file hybridpressionmodulempi.h.

Constructor & Destructor Documentation

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 = std::cout 
)

Definition at line 15 of file hybridpressionmodulempi.cpp.

00016   :HybridPressionModule(mesh,fPrescribedVelocities,fPrescribedPression,K,fMobT,fGravSource,debugLevel,false),
00017    m_out(out)
00018 {
00019   assert(m_theta >= 0.0 && m_theta <=1.0);
00020   m_max_tol = maxTol;
00021   m_meshOverlap = meshOverlap;
00022   setCommunication();
00023   m_xSize = Xsize;
00024   m_b=b;
00025   m_sol=0.0;
00026   m_solPrev=m_sol;
00027   m_theta=theta;
00028   if (meshOverlap<=0)
00029     throw new Exception("Invalid meshOverlap parameter");
00030 
00031 #ifdef DEBUG
00032   checkBC(fPrescribedVelocities,fPrescribedPression);
00033 #endif
00034 
00035   //Retrieve this code in the final release
00036     VecDouble Vv(m_mesh.numVertices());
00037     HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 
00038 
00039     
00040 }

HybridPressionModuleMPI::~HybridPressionModuleMPI (  )  [virtual]

Definition at line 271 of file hybridpressionmodulempi.cpp.

00272 {
00273 
00274 }

Member Function Documentation

void HybridPressionModuleMPI::addRobinBC ( const VecDouble &  pB,
const VecIndex &  fI,
const VecDouble &  vK,
const VecDouble &  S1,
Function1D &  fMobT,
SparseMatrix< double > &  A 
) [protected]

Definition at line 404 of file hybridpressionmodulempi.cpp.

00405 {
00406   assert(vB.size() == fI.size());
00407   assert(vK.size() == m_mesh.numCells());
00408   assert(S1.size() == m_mesh.numCells());
00409   assert(A.m() == A.n() && A.m() == m_mesh.numCells());
00410 
00411   
00412   for (unsigned i=0;i<vB.size();i++)
00413   {
00414     OrthoMesh::Face_It face = m_mesh.get_face(fI[i]);
00415     assert(face->at_boundary());
00416     assert(face->getNormalOrientation() == OrthoMesh::NORMAL_X);
00417     unsigned c1,c2;
00418     face->getAdjCellIndices(c1,c2);
00419     double fArea = face->area(); 
00420     if (c2 != OrthoMesh::INVALID_INDEX) //left boundary
00421     {
00422       double K = vK(c2)*fMobT(S1(c2));
00423       double BC =  2*K*fArea/(m_mesh.get_dx() + 2*K*vB(i));
00424       A.add(c2,c2,BC);
00425     }
00426     if (c1 != OrthoMesh::INVALID_INDEX)//right boundary
00427     {
00428       double K = vK(c1)*fMobT(S1(c1));
00429       double BC = 2*K*fArea/(m_mesh.get_dx() + 2*K*vB(i));
00430       A.add(c1,c1,BC);
00431     }
00432     
00433     
00434   }
00435   
00436 }

void HybridPressionModuleMPI::addRobinBC ( const VecDouble &  v1,
const VecDouble &  pB,
const VecIndex &  fI,
const VecDouble &  vK,
const VecDouble &  S1,
Function1D &  fMobT,
VecDouble &  vRHS 
) [protected]

Definition at line 366 of file hybridpressionmodulempi.cpp.

00367 {
00368   assert(vB.size() == fI.size());
00369   assert(v1.size() == fI.size());
00370   assert(vK.size() == m_mesh.numCells());
00371   assert(S1.size() == m_mesh.numCells());
00372   assert(vRHS.size() == m_mesh.numCells());
00373  
00374   for (unsigned i=0;i<vB.size();i++)
00375   {
00376     OrthoMesh::Face_It face = m_mesh.get_face(fI[i]);
00377     if (!face->at_boundary())
00378     {
00379       std::cout << "throw exception\n" << std::endl;
00380       throw new Exception("Problem at face %d",fI[i]);
00381 
00382     }
00383     assert(face->at_boundary());
00384     assert(face->getNormalOrientation() == OrthoMesh::NORMAL_X);
00385     unsigned c1,c2;
00386     face->getAdjCellIndices(c1,c2);
00387     double fArea = face->area(); 
00388     if (c2 != OrthoMesh::INVALID_INDEX)
00389     {
00390       double K = vK(c2)*fMobT(S1(c2));
00391       vRHS(c2)+= K*2*v1(i)*fArea/(m_mesh.get_dx() + 2*K*vB(i));
00392     }
00393     if (c1 != OrthoMesh::INVALID_INDEX)
00394     {
00395       double K = vK(c1)*fMobT(S1(c1));
00396       vRHS(c1)+= K*2*v1(i)*fArea/(m_mesh.get_dx() + 2*K*vB(i));
00397     }
00398     
00399     
00400   }
00401 }

void HybridPressionModuleMPI::assemblySystem ( TransportBase &  trans,
SparseMatrix< double > &  A,
VecDouble &  B,
VecDouble &  K,
const VecDouble &  S1,
Function1D &  fMobT,
FunctionOfCellFields &  fMobGrav 
) [protected]

Definition at line 536 of file hybridpressionmodulempi.cpp.

00537 {
00538   //SET AND GET  MPI VARIABLES 
00539   MPI_Request req[4];
00540   MPI_Status stats[4];
00541   req[0]=req[1]=req[2]=req[3]=MPI_REQUEST_NULL;
00542    
00543   int rank = NetMPI::rank();
00544   int next = rank+1;
00545   int prev = rank-1;
00546 
00547   
00548   if (rank%2 == 0)
00549   {
00550 
00551     //First of all we need to get the B parameter used in the domain decomposition
00552 
00553     //Send B to left and right domain
00554     getBParameter(BL_Snd,K,S1,fMobT,fLeftBCSndMap);
00555     NetMPI::Isend(BL_Snd,prev,B_PARAMETER,req[0]);
00556     getBParameter(BR_Snd,K,S1,fMobT,fRightBCSndMap);
00557     NetMPI::Isend(BR_Snd,next,B_PARAMETER,req[1]);
00558 
00559     //Receive B from left  and right  domain
00560     NetMPI::Irecv(BL_Rcv,prev,B_PARAMETER,req[2]);
00561     NetMPI::Irecv(BR_Rcv,next,B_PARAMETER,req[3]);
00562 
00563     //Set the gravity mobility function defined in the cells
00564     //with the solutions of the transport module
00565     for (unsigned fieldId=0;fieldId<fMobGrav.numFields();fieldId++) //Set the fields of the gravity source term
00566       fMobGrav.setField(trans.getSolutionAtCells(fieldId),fieldId);
00567 
00568     //assembly system
00569     HybridPressionModule::assemblySystem(A,B,K,S1,fMobT,fMobGrav);
00570 
00571 
00572     //Sync point reached
00573     //after this point all processors have the B parameter from its neighbors
00574     MPI_Waitall(4,req,stats);
00575 
00576     addRobinBC(BL_Rcv,fLeftBCRcvMap,K,S1,fMobT,A);
00577     addRobinBC(BR_Rcv,fRightBCRcvMap,K,S1,fMobT,A);
00578 
00579   
00580   }
00581   else if (rank%2 == 1)
00582   {
00583 
00584     //First of all we need to get the B parameter used in the domain decomposition
00585 
00586     //Receive B from right and left domain
00587     NetMPI::Irecv(BR_Rcv,next,B_PARAMETER,req[0]);
00588     NetMPI::Irecv(BL_Rcv,prev,B_PARAMETER,req[1]);
00589 
00590     //Send B to right and left domain
00591     getBParameter(BR_Snd,K,S1,fMobT,fRightBCSndMap);
00592     NetMPI::Isend(BR_Snd,next,B_PARAMETER,req[3]);
00593     getBParameter(BL_Snd,K,S1,fMobT,fLeftBCSndMap);
00594     NetMPI::Isend(BL_Snd,prev,B_PARAMETER,req[2]);
00595 
00596     //Set the gravity mobility function defined in the cells
00597     //with the solutions of the transport module
00598     for (unsigned fieldId=0;fieldId<fMobGrav.numFields();fieldId++) //Set the fields of the gravity source term
00599       fMobGrav.setField(trans.getSolutionAtCells(fieldId),fieldId);
00600 
00601     //assembly system
00602     HybridPressionModule::assemblySystem(A,B,K,S1,fMobT,fMobGrav);
00603 
00604 
00605     //Sync point reached
00606     //after this point all processors have the B parameter from its neighbors
00607     MPI_Waitall(4,req,stats);
00608 
00609     //With the B parameter we can add the robin conditions to the system matrix.
00610     addRobinBC(BL_Rcv,fLeftBCRcvMap,K,S1,fMobT,A);
00611     addRobinBC(BR_Rcv,fRightBCRcvMap,K,S1,fMobT,A);
00612   }
00613 }

void HybridPressionModuleMPI::checkBC ( Function3D &  fPrescribedVelocities,
Function3D &  fPrescribedPression 
) [protected]

Check if all boundary faces have boundary conditions

Definition at line 88 of file hybridpressionmodulempi.cpp.

00089 {
00090   Point3D p;  
00091 
00092   OrthoMesh::Cell_It cell = m_mesh.begin_cell(),
00093     endc = m_mesh.end_cell();
00094 
00095   
00096    //For eahc cell in the Mesh
00097    for(;cell != endc;cell++)
00098    {
00099      //For each boundary face of the cell
00100      for (unsigned int face_no = 0; face_no < GeometryInfo<DIM>::faces_per_cell;++face_no)
00101      {
00102        FaceDirection3D fDir = (FaceDirection3D) face_no;
00103        OrthoMesh::Face_It face = cell->face(fDir);
00104        if (face->at_boundary())
00105        {
00106          bool found=false;
00107          face->barycenter(p);
00108          if (fPrescribedPression.isInDomain(p))
00109          {
00110            found=true;
00111          }
00112          if (fPrescribedVelocities.isInDomain(p))
00113          {
00114            if (found)
00115              throw new Exception("Boundary conditions clash");
00116            else
00117              found=true;
00118          }
00119          if ( (std::find(fLeftBCRcvMap.begin(),fLeftBCRcvMap.end(),face->index()) !=fLeftBCRcvMap.end()) ||
00120               (std::find(fRightBCRcvMap.begin(),fRightBCRcvMap.end(),face->index()) !=fRightBCRcvMap.end()))
00121          {
00122            if (found)
00123              throw new Exception("Boundary conditions clash");
00124            else
00125              found=true;
00126          }
00127          if (!found)
00128          {
00129            Point3D p;
00130            face->barycenter(p); 
00131            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]);
00132          }
00133            
00134        }
00135      }
00136    }
00137 }

int HybridPressionModuleMPI::checkConvergence ( double  tol  )  [protected]

Definition at line 615 of file hybridpressionmodulempi.cpp.

00616 {
00617   int bContinue;
00618   HDF5OrthoWriter::getHDF5OrthoWriter().setVariable("ItCnt",static_cast<double>(m_itCnt));
00619 
00620   if (NetMPI::rank() != 0)
00621   {
00622     MPI_Gather(&tol,1,MPI_DOUBLE,NULL,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
00623     MPI_Bcast(&bContinue,1,MPI_INT,0,MPI_COMM_WORLD);
00624   }
00625   else
00626   {
00627     m_vTol(0)=tol;
00628     MPI_Gather(&tol,1,MPI_DOUBLE,&(m_vTol(0)),1,MPI_DOUBLE,0,MPI_COMM_WORLD);
00629     bContinue=false;
00630     printf("%d Tolerance: ",m_itCnt++);
00631     for (unsigned i=0;i<m_vTol.size();i++)
00632     {
00633       printf("%g ",m_vTol(i));
00634       if (m_vTol(i) >= m_max_tol)
00635       {
00636         bContinue = true;
00637         break;
00638       }
00639     }
00640     printf("\n");
00641     MPI_Bcast(&bContinue,1,MPI_INT,0,MPI_COMM_WORLD);
00642   }
00643   return bContinue;
00644 }

void HybridPressionModuleMPI::exchangeData ( VecDouble &  lSnd,
VecDouble &  lRcv,
VecDouble &  rSnd,
VecDouble &  rRcv 
) [protected]

Exchange data between neighbor domains

Parameters:
lSnd Data to send to left domain
lRcv Data to be received from left domain
rSnd Data to be sent to the right domain
rRcv Data to be received from right domain
Returns:

Definition at line 283 of file hybridpressionmodulempi.cpp.

00284 {
00285 
00286 
00287   int rank = NetMPI::rank();
00288   int next=rank+1;
00289   int prev=rank-1;
00290   MPI_Request reqs[4];
00291   reqs[0]=reqs[1]=reqs[2]=reqs[3]=MPI_REQUEST_NULL;
00292   
00293   MPI_Status stats;
00294 
00295   assert(!(rank == 0 && lSnd.size() != 0));
00296   assert(!(NetMPI::isLastProcess() && rSnd.size() != 0));
00297   
00298   if (rank%2 == 0) //black domain
00299   {
00300     if (rank!= 0)
00301       MPI_Send(&(lSnd(0)),lSnd.size(),MPI_DOUBLE,prev,BLACK_TO_RED,MPI_COMM_WORLD);
00302     if (!NetMPI::isLastProcess())
00303       MPI_Send(&(rSnd(0)),rSnd.size(),MPI_DOUBLE,next,BLACK_TO_RED,MPI_COMM_WORLD);
00304 
00305     if (rank!=0)
00306       MPI_Recv(&(lRcv(0)),lRcv.size(),MPI_DOUBLE,prev,RED_TO_BLACK,MPI_COMM_WORLD,&stats);
00307     if (!NetMPI::isLastProcess())
00308       MPI_Recv(&(rRcv(0)),rRcv.size(),MPI_DOUBLE,next,RED_TO_BLACK,MPI_COMM_WORLD,&stats);
00309   }
00310   else
00311   {
00312     if (!NetMPI::isLastProcess())
00313       MPI_Recv(&(rRcv(0)),rRcv.size(),MPI_DOUBLE,next,BLACK_TO_RED,MPI_COMM_WORLD,&stats);
00314     MPI_Recv(&(lRcv(0)),lRcv.size(),MPI_DOUBLE,prev,BLACK_TO_RED,MPI_COMM_WORLD,&stats);
00315 
00316     if (!NetMPI::isLastProcess())
00317       MPI_Send(&(rSnd(0)),rSnd.size(),MPI_DOUBLE,next,BLACK_TO_RED,MPI_COMM_WORLD);
00318     MPI_Send(&(lSnd(0)),lSnd.size(),MPI_DOUBLE,prev,BLACK_TO_RED,MPI_COMM_WORLD);
00319 
00320     
00321 
00322   }
00323 
00324 }

void HybridPressionModuleMPI::getBParameter ( VecDouble &  vB,
const VecDouble &  K,
const VecDouble &  S1,
Function1D &  fMobT,
const VecIndex &  fI 
) [protected]

Definition at line 327 of file hybridpressionmodulempi.cpp.

00328 {
00329   assert(vB.size() == fI.size());
00330   for (unsigned i=0;i<vB.size();i++)
00331   {
00332     OrthoMesh::Face_It face = m_mesh.get_face(fI[i]);
00333     assert(!face->at_boundary());
00334     unsigned c1,c2;
00335     face->getAdjCellIndices(c1,c2);
00336 
00337     double K1 = K(c1)*fMobT(S1(c1));
00338     double K2 = K(c2)*fMobT(S1(c2));
00339      
00340     double Keff = 2*K1*K2/(K1 + K2);
00341     vB(i)=m_b*m_xSize/Keff;
00342   }
00343 
00344 }

void HybridPressionModuleMPI::getNormalVelocitiesFromPressure ( VecDouble &  Vq,
const VecDouble &  K,
const VecDouble &  vP,
const VecDouble &  S1,
Function1D &  fMobT,
FunctionOfCellFields &  fMobGrav,
const VecDouble &  LRobinBC,
const VecDouble &  RRobinBC 
) [protected]

Definition at line 502 of file hybridpressionmodulempi.cpp.

00503 {
00504   HybridPressionModule::getNormalVelocitiesFromPressure(Vq,K,vP,S1,fMobT,fMobGrav);
00505   assert(RRobinBC.size() ==0 || RRobinBC.size() == fRightBCRcvMap.size());
00506   assert(LRobinBC.size() ==0 || LRobinBC.size() == fLeftBCRcvMap.size());
00507 
00508   for (unsigned i=0;i<LRobinBC.size();i++)
00509   {
00510     int faceIndex = fLeftBCRcvMap[i];
00511     OrthoMesh::Face_It face = m_mesh.get_face(faceIndex);
00512     unsigned  c2 = face->getNegCell();
00513     assert(face->at_boundary());
00514     assert(face->getNormalOrientation() == OrthoMesh::NORMAL_X);
00515     assert(c2!=OrthoMesh::INVALID_INDEX);
00516     double k = K(c2)*fMobT(S1(c2));
00517     double p = vP(c2);
00518     Vq(faceIndex) = - 2*k*(p-LRobinBC(i))/(m_mesh.get_dx() + 2*k*BL_Rcv(i));
00519   }
00520   
00521   for (unsigned i=0;i<RRobinBC.size();i++)
00522   {
00523     int faceIndex = fRightBCRcvMap[i];
00524     OrthoMesh::Face_It face = m_mesh.get_face(faceIndex);
00525     unsigned c1 = face->getPosCell();
00526     assert(face->at_boundary());
00527     assert(face->getNormalOrientation() == OrthoMesh::NORMAL_X);
00528     assert(c1!=OrthoMesh::INVALID_INDEX);
00529     double k = K(c1)*fMobT(S1(c1));
00530     double p = vP(c1);
00531     Vq(faceIndex) = - 2*k*(RRobinBC(i)-p)/(m_mesh.get_dx() + 2*k*BR_Rcv(i));
00532   }
00533 }

void HybridPressionModuleMPI::getRobinBC ( VecDouble &  vRobin,
const VecDouble &  K,
const VecDouble &  vP,
const VecDouble &  S1,
Function1D &  fMobT,
const VecIndex &  fI,
bool  bNegateB 
) [protected]

Definition at line 347 of file hybridpressionmodulempi.cpp.

00348 {
00349   assert(vRobin.size() == fI.size());
00350   double b = bNegateB ? -m_b*m_xSize : m_b*m_xSize;
00351   for (unsigned i=0;i<vRobin.size();i++)
00352   {
00353     OrthoMesh::Face_It face = m_mesh.get_face(fI[i]);
00354     assert(!face->at_boundary());
00355     unsigned c1,c2;
00356     face->getAdjCellIndices(c1,c2);
00357     double p1 = vP(c1);
00358     double p2 = vP(c2);
00359     double K1 = K(c1)*fMobT(S1(c1));
00360     double K2 = K(c2)*fMobT(S1(c2));
00361     vRobin(i) = +b*(p2-p1)*face->areaPerCellVol() + (K1*p1 + K2*p2)/(K1 + K2);
00362   }
00363 }

virtual void HybridPressionModuleMPI::initIteration (  )  [protected, pure virtual]

Implemented in HybridPressionModuleAGM_MPI, and HybridPressionModuleUMFPACK_MPI.

void HybridPressionModuleMPI::iterate ( TransportBase &  trans  )  [virtual]

Definition at line 439 of file hybridpressionmodulempi.cpp.

00440 {
00441   int rank = NetMPI::rank();
00442 
00443   if (rank%2 == 0)
00444   {
00445     iterateRed(trans);
00446 
00447   }
00448   else if (rank%2 == 1)
00449   {
00450     iterateBlack(trans);
00451   }
00452     
00453 }

void HybridPressionModuleMPI::iterateBlack ( TransportBase &  trans  )  [protected, virtual]

The iteration of black nodes. The black nodes are the process with odd indices |R|B|R|B|......

Parameters:
trans 
Returns:

Definition at line 206 of file hybridpressionmodulempi.cpp.

00207 {
00208   //SET AND GET  MPI VARIABLES 
00209   MPI_Request _req[4];
00210   MPI_Status stats[4];
00211   _req[0]=_req[1]=_req[2]=_req[3]=MPI_REQUEST_NULL;
00212 
00213   double tol;
00214   int rank = NetMPI::rank();
00215   int next = rank+1;
00216   int prev = rank-1;
00217 
00218   const VecDouble &S1 = trans.getSolutionAtCells(0);
00219   assemblySystem(trans,m_A,m_B,m_K,S1,m_fMobT,m_fMobGrav);
00220   initIteration();
00221   do {
00222 
00223     m_RHS = m_B;
00224     addRobinBC(LBC_Rcv,BL_Rcv,fLeftBCRcvMap,m_K,S1,m_fMobT,m_RHS);
00225     addRobinBC(RBC_Rcv,BR_Rcv,fRightBCRcvMap,m_K,S1,m_fMobT,m_RHS);
00226     m_solPrev=m_sol;
00227     solve();
00228 
00229     //Receive right and left  robin condition
00230     RBC_Prev.swap(RBC_Rcv);
00231     LBC_Prev.swap(LBC_Rcv);
00232     NetMPI::Irecv(RBC_Rcv,next,ROBIN_BC_MESSAGE,_req[0]);
00233     NetMPI::Irecv(LBC_Rcv,prev,ROBIN_BC_MESSAGE,_req[1]);
00234 
00235     //Send right and left robin condition
00236     getRobinBC(RBC_Snd,m_K,m_sol,S1,m_fMobT,fRightBCSndMap,false);
00237     NetMPI::Isend(RBC_Snd,next,ROBIN_BC_MESSAGE,_req[3]);
00238     getRobinBC(LBC_Snd,m_K,m_sol,S1,m_fMobT,fLeftBCSndMap,true);
00239     NetMPI::Isend(LBC_Snd,prev,ROBIN_BC_MESSAGE,_req[4]);
00240     MPI_Waitall(4,_req,stats);
00241 
00242     
00243     if (rank != 0)
00244       LBC_Rcv.sadd((1.0-m_theta),m_theta,LBC_Prev);
00245     if (!NetMPI::isLastProcess())
00246       RBC_Rcv.sadd((1.0-m_theta),m_theta,RBC_Prev);
00247 
00248     //tol = std::max(NumericMethods::relError(LBC_Rcv,LBC_Prev),
00249     //             NumericMethods::relError(RBC_Rcv,RBC_Prev));
00250 
00251     tol = NumericMethods::relError(m_solPrev,m_sol);
00252 
00253 
00254   }while (checkConvergence(tol));
00255 
00256   getNormalVelocitiesFromPressure(m_Vq,m_K,m_sol,S1,m_fMobT,m_fMobGrav,LBC_Rcv,RBC_Rcv);
00257   
00258 }

void HybridPressionModuleMPI::iterateRed ( TransportBase &  trans  )  [protected, virtual]

The iteration of red nodes. The red nodes are the process with even indices counting from 0 |R|B|R|B|......

Parameters:
trans 
Returns:

Definition at line 145 of file hybridpressionmodulempi.cpp.

00146 {
00147   //SET AND GET  MPI VARIABLES 
00148   MPI_Request _req[4];
00149   MPI_Status stats[4];
00150   _req[0]=_req[1]=_req[2]=_req[3]=MPI_REQUEST_NULL;
00151 
00152   double tol;
00153   int rank = NetMPI::rank();
00154   int next = rank+1;
00155   int prev = rank-1;
00156 
00157 
00158   const VecDouble &S1 = trans.getSolutionAtCells(0);
00159 
00160   assemblySystem(trans,m_A,m_B,m_K,S1,m_fMobT,m_fMobGrav);
00161   initIteration();
00162   
00163 
00164   m_itCnt=0;
00165   do {
00166     m_RHS = m_B;
00167     addRobinBC(LBC_Rcv,BL_Rcv,fLeftBCRcvMap,m_K,S1,m_fMobT,m_RHS);
00168     addRobinBC(RBC_Rcv,BR_Rcv,fRightBCRcvMap,m_K,S1,m_fMobT,m_RHS);
00169     m_solPrev=m_sol;
00170     solve();
00171 
00172     //Send left and right robin condition
00173     getRobinBC(LBC_Snd,m_K,m_sol,S1,m_fMobT,fLeftBCSndMap,true);
00174     NetMPI::Isend(LBC_Snd,prev,ROBIN_BC_MESSAGE,_req[0]);
00175     getRobinBC(RBC_Snd,m_K,m_sol,S1,m_fMobT,fRightBCSndMap,false);
00176     NetMPI::Isend(RBC_Snd,next,ROBIN_BC_MESSAGE,_req[1]);
00177 
00178     //Receive left and right  robin condition
00179     RBC_Prev.swap(RBC_Rcv);
00180     LBC_Prev.swap(LBC_Rcv);
00181     NetMPI::Irecv(LBC_Rcv,prev,ROBIN_BC_MESSAGE,_req[2]);
00182     NetMPI::Irecv(RBC_Rcv,next,ROBIN_BC_MESSAGE,_req[3]);
00183     MPI_Waitall(4,_req,stats);
00184  
00185     if (rank != 0)
00186       LBC_Rcv.sadd((1.0-m_theta),m_theta,LBC_Prev);
00187     if (!NetMPI::isLastProcess())
00188       RBC_Rcv.sadd((1.0-m_theta),m_theta,RBC_Prev);
00189 
00190     //    tol = std::max(NumericMethods::relError(LBC_Rcv,LBC_Prev),
00191     //             NumericMethods::relError(RBC_Rcv,RBC_Prev));
00192     tol = NumericMethods::relError(m_solPrev,m_sol);
00193   } while (checkConvergence(tol));
00194 
00195   getNormalVelocitiesFromPressure(m_Vq,m_K,m_sol,S1,m_fMobT,m_fMobGrav,LBC_Prev,RBC_Prev);
00196 
00197 }

void HybridPressionModuleMPI::printOutput (  )  [virtual]

Definition at line 458 of file hybridpressionmodulempi.cpp.

00459 {
00460   if (_debugLevel > 0)
00461   {
00462     VecDouble v(m_mesh.numCells());
00463     VecDouble Vv(m_mesh.numVertices());
00464     HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 
00465 
00466     this->getPressionAtCells(v);
00467     m_mesh.projectCentralValuesAtVertices(v,Vv);
00468     hdf5.writeScalarField(Vv,"P");
00469 
00470     //    GnuPlotAnim &plot = GnuPlotAnim::getGnuPlotAnim();
00471     //plot.plotVerticeValuesAtFixedZ(m_mesh.getTriangulation(),Vv,"pression", "with lines",250);
00472     //plot.NextScene();
00473 
00474   
00475     if (_debugLevel > 1)
00476     {
00477       Matrix Mvel(m_mesh.numCells(),3);
00478       this->getVelocitiesAtCells(Mvel);
00479       if (_debugLevel > 2)
00480       {
00481         printf("\nVelocities at cells' barycenters\n\n");
00482         Mvel.print_formatted(m_out,3,true,0,"0");
00483       }
00484       VecDouble v(m_mesh.numCells());
00485       for (unsigned i=0;i<v.size();i++)
00486       {
00487         v(i) = Mvel(i,0); 
00488       }
00489       VecDouble minVel(3);
00490       VecDouble maxVel(3);
00491       NumericMethods::getMinMaxValueByColumn(Mvel,minVel,maxVel);
00492       printf("Min Vel <%g,%g,%g> - Max Vel <%g,%g,%g>\n",minVel(0),minVel(1),minVel(2),maxVel(0),maxVel(1),maxVel(2));
00493   
00494       m_mesh.projectCentralValuesAtVertices(v,Vv);
00495       hdf5.writeScalarField(Vv,"Vel");
00496     }
00497   }
00498 
00499 }

void HybridPressionModuleMPI::setCommunication (  )  [protected]

Definition at line 42 of file hybridpressionmodulempi.cpp.

00043 {
00044   if (NetMPI::rank()!=0)
00045   {
00046     m_mesh.getFacesInYZSlab(m_mesh.getP()[0],fLeftBCRcvMap);
00047     m_mesh.getFacesInYZSlab(m_mesh.getP()[0]+2*m_meshOverlap*m_mesh.getDX()[0],fLeftBCSndMap);
00048 
00049     //get B parameter to send to the left domain
00050     BL_Snd.reinit(fLeftBCRcvMap.size());
00051     BL_Rcv.reinit(fLeftBCRcvMap.size());
00052     LBC_Snd.reinit(fLeftBCRcvMap.size());
00053     LBC_Rcv.reinit(fLeftBCRcvMap.size());
00054     LBC_Prev.reinit(fLeftBCRcvMap.size());
00055 
00056     LBC_Rcv=0.0;
00057     LBC_Prev=0.0;
00058     
00059   }
00060   if (!NetMPI::isLastProcess())
00061   {
00062     m_mesh.getFacesInYZSlab(m_mesh.getQ()[0],fRightBCRcvMap);
00063     m_mesh.getFacesInYZSlab(m_mesh.getQ()[0]-2*m_meshOverlap*m_mesh.getDX()[0],fRightBCSndMap);
00064 
00065     //get B parameter to send to the right domain
00066     BR_Snd.reinit(fRightBCRcvMap.size());
00067     BR_Rcv.reinit(fRightBCRcvMap.size());
00068     RBC_Snd.reinit(fRightBCRcvMap.size());
00069     RBC_Rcv.reinit(fRightBCRcvMap.size());
00070     RBC_Prev.reinit(fRightBCRcvMap.size());
00071     RBC_Rcv=0.0;
00072     RBC_Prev=0.0;
00073   }
00074   if (NetMPI::rank() == 0)
00075   {
00076     m_vTol.reinit(NetMPI::nProcess());
00077   }
00078   
00079   
00080 
00081 }

virtual void HybridPressionModuleMPI::solve (  )  [protected, pure virtual]

Implemented in HybridPressionModuleAGM_MPI, and HybridPressionModuleUMFPACK_MPI.

Member Data Documentation

VecIndex HybridPressionModuleMPI::fLeftBCRcvMap [protected]

The sequence of faces indices containing robin boundary conditions to be received from the left subdomain

Definition at line 22 of file hybridpressionmodulempi.h.

VecIndex HybridPressionModuleMPI::fLeftBCSndMap [protected]

The sequence of faces indices containing robin boundary conditions to be sent to the left subdomain

Definition at line 24 of file hybridpressionmodulempi.h.

VecIndex HybridPressionModuleMPI::fRightBCRcvMap [protected]

The sequence of faces indices containing robin boundary conditions to be received from the right subdomain

Definition at line 23 of file hybridpressionmodulempi.h.

VecIndex HybridPressionModuleMPI::fRightBCSndMap [protected]

The sequence of faces indices containing robin boundary conditions to be sent to the right subdomain

Definition at line 25 of file hybridpressionmodulempi.h.

double HybridPressionModuleMPI::m_b [protected]

Parameter that is part of the Beta parameter definition

Definition at line 30 of file hybridpressionmodulempi.h.

unsigned HybridPressionModuleMPI::m_itCnt [protected]

Definition at line 32 of file hybridpressionmodulempi.h.

double HybridPressionModuleMPI::m_max_tol [protected]

Max error of the convergence

Definition at line 50 of file hybridpressionmodulempi.h.

int HybridPressionModuleMPI::m_meshOverlap [protected]

How many elements along the X plane are shared by adjacent subdomains

Definition at line 28 of file hybridpressionmodulempi.h.

std::ostream& HybridPressionModuleMPI::m_out [protected]

Definition at line 21 of file hybridpressionmodulempi.h.

VecDouble HybridPressionModuleMPI::m_RHS [protected]

Definition at line 33 of file hybridpressionmodulempi.h.

VecDouble HybridPressionModuleMPI::m_solPrev [protected]

Definition at line 34 of file hybridpressionmodulempi.h.

double HybridPressionModuleMPI::m_theta [protected]

Parameter Theta for the interaction

Definition at line 31 of file hybridpressionmodulempi.h.

VecDouble HybridPressionModuleMPI::m_vTol [protected]

Vector to contain the errors of all processes. Relevant only for the root process

Definition at line 49 of file hybridpressionmodulempi.h.

double HybridPressionModuleMPI::m_xSize [protected]

X size of the domain

Definition at line 29 of file hybridpressionmodulempi.h.

The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Sun Apr 8 23:13:12 2012 for CO2INJECTION by  doxygen 1.6.3