parallel/pico/DomainCoupler.cpp

// -*- C++ -*-
//
//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
//                                   Fehmi Cirak
//                        California Institute of Technology
//                           (C) 2004 All Rights Reserved
//
//  <LicenseText>
//
//  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// EDITED BY MATT GALLOWAY FOR USE WITH DAPTA

#include "DAPTA.hpp"

#include <limits>
#include <cassert>
#include <cstdlib>

namespace DAPTA { // Define namespace DAPTA

namespace pico { // Define namespace pico

template <typename V>
void DomainCoupler<V>::registerSubDomain(const SolvType& type, const CoordArrayType lo, 
                                 const CoordArrayType up, const double& maxElemSize) {
   VertexType loV(lo);
   VertexType upV(up);
   
   SubDomain<V> *sd = new SubDomain<V>(loV, upV, maxElemSize);
   
   int rank;
   MPI_Comm_rank(_pComm, &rank);
   
   _subDomains[type][rank] = sd;
}

template <typename V>
void DomainCoupler<V>::registerSubDomain(const SolvType& type, const CoordType * const coor, 
                               const int& nodes, const double& maxElemSize) {
   std::numeric_limits<CoordType> doubleLimits;
   CoordType max = doubleLimits.max();
   CoordType min = doubleLimits.min();
   
   CoordArrayType lo;
   CoordArrayType up;
   
   for(int i=0; i<V::dimension; ++i) {
      lo[i] = max;
      up[i] = min;
   }
   
   for(int i=0; i<nodes; ++i) {
      for(int j=0; j<V::dimension; ++j) {
         lo[j] = std::min(lo[j], coor[V::dimension*i+j]);
         up[j] = std::max(up[j], coor[V::dimension*i+j]);
      }
   }
   
   registerSubDomain(type, lo, up, maxElemSize);
}

template <typename V>
void DomainCoupler<V>::gatherSubDomains() {
   const int numVarHalf = 8;
   const int numVar = 2*numVarHalf;
   double myRBuf[numVar];
   
   int rank;
   MPI_Comm_rank(_pComm, &rank);    
   
   VertexType loV, upV;
   for (int i=0; i<2; ++i) {
      
      _SubDomainIt its, itse = _subDomains[i].end();
      its = _subDomains[i].find(rank);
      
      int start = i*numVarHalf;
      
      if (its==itse) 
         myRBuf[start++] = -1;
      else 
         myRBuf[start++] = i;
      
      if (its==itse) continue;
      
      SubDomain<V> *sd = its->second;
      
      loV = sd->lowerPoint();
      upV = sd->upperPoint();
      
      for(int j=0; j<V::dimension; j++) {
         myRBuf[start++] = loV.coords[j];
         myRBuf[start++] = upV.coords[j];
      }
      myRBuf[start++] = sd->maxElemSize();
   }
   
   int numDomains;
   MPI_Comm_size(_pComm, &numDomains);
   
   double *rBuf = new double[numDomains*numVar];
   
   MPI_Allgather(myRBuf, numVar, MPI_DOUBLE, rBuf, numVar, MPI_DOUBLE, _pComm); 
   
   // create the subdomains on other processors
   for (int j=0; j<numDomains; ++j) {
      
      if (j==rank) continue;
      
      for (int i=0; i<2; ++i) {
         int start = j*numVar+i*numVarHalf;
         
         int domainType = static_cast<int>(rBuf[start++]);
         if (domainType == -1 ) continue;
         
         for(int k=0; k<V::dimension; k++) {
            loV.coords[k] = rBuf[start++];
            upV.coords[k] = rBuf[start++];
         }
         
         double elemSize = rBuf[start++];
         
         SubDomain<V> *sd = new SubDomain<V>(loV, upV, elemSize);
         
         _subDomains[domainType][j] = sd;
      }
   }
   
   if (rBuf) delete[] rBuf;
}

template <typename V>
void DomainCoupler<V>::prepareCommunicationsInter() {
   // gather the bounding box data of all subdomains
   gatherSubDomains();
   
   int rank, otherType[2]={1, 0};
   MPI_Comm_rank(_pComm, &rank);  
   
   // check if there are two types of subDomains
   _SubDomainIt itt, itb[2], ite[2];
   for (int type=0; type<2; ++type) {
      itb[type] = _subDomains[type].begin();
      ite[type] = _subDomains[type].end();
      if (itb[type]==ite[type]) return;
   }
   
   // find all relevant subdomains of opposite type
   for (int type=0; type<2; ++type) {  
      
      itt = _subDomains[type].find(rank);
      if (itt==ite[type]) continue;
      SubDomain<V>* myDomain = itt->second;
      
      _SubDomainIt its = itb[otherType[type]], itse = ite[otherType[type]];
      
      // make my bounding box larger by the requested element overlapsize
      // think about it
      
      for ( ; its!=itse; ++its) {       
         SubDomain<V> *sd = its->second;      
         if (sd->checkOverlap(myDomain)) {
            _relevSubDomainsInter[otherType[type]][its->first] = sd;
         }
      }  
   }
}

template <typename V>
void DomainCoupler<V>::prepareCommunicationsIntra(const SolvType& type) {
   // gather the bounding box data of all subdomains
   gatherSubDomains();
   
   int rank;
   MPI_Comm_rank(_pComm, &rank);  
   
   // check if there is a solver of the type type
   _SubDomainIt it, ite=_subDomains[type].end();
   
   it = _subDomains[type].find(rank);
   if (it==ite) return;
   SubDomain<V> *myDomain = it->second;
   
   // make my bounding box larger by the requested element overlapsize
   // this might be later neccessary for shells
   
   it = _subDomains[type].begin();
   for ( ; it!=ite; ++it) {
      SubDomain<V> *sd = it->second;
      if (sd==myDomain) continue;
      
      if (sd->checkOverlap(myDomain)) {
         _relevSubDomainsIntra[type][it->first] = sd;
      }
   }
}

template <typename V>
void DomainCoupler<V>::exchangeBufferSizesIntra(const SolvType& type) {
   int rank;        
   MPI_Comm_rank(_pComm, &rank);  
   
   _SubDomainIt its, itse = _subDomains[type].end();
   its = _subDomains[type].find(rank);
   
   if (its==itse) return;    
   
   SubDomain<V>* myDomain = its->second;
   
   int mySBuf = myDomain->sendMsgSize();
   
   _SubDomainIt it, itb = _relevSubDomainsIntra[type].begin(), ite = _relevSubDomainsIntra[type].end();
   
   MPI_Status status;
   int tag=0, rBuf=0, myRBuf;
   for (it=itb; it!=ite; ++it) {
      MPI_Sendrecv(&mySBuf, 1, MPI_INT, it->first, tag, &myRBuf, 1, MPI_INT, it->first, tag, _pComm, &status); 
      (it->second)->setSendMsgSize(myRBuf);
      rBuf += myRBuf;
   }
   
   _subDomains[type][rank]->allocateReceiveBuffer(rBuf);     
   
   assert(rBuf>0);
}

template <typename V>
void DomainCoupler<V>::exchangeSubdomainDataIntra(const SolvType& type) {    
   exchangeBufferSizesIntra(type);
   
   int rank;
   MPI_Comm_rank(_pComm, &rank);     
   SubDomain<V>* myDomain = _subDomains[type][rank];    
   char *recvBuf = myDomain->recvBuffer();
   
   MPI_Request *request = 
      (MPI_Request *) std::malloc(_relevSubDomainsIntra[type].size()*sizeof(MPI_Request));
   MPI_Status *status = 
      (MPI_Status *) std::malloc(_relevSubDomainsIntra[type].size()*sizeof(MPI_Status));
   
   _SubDomainIt it, itb = _relevSubDomainsIntra[type].begin();
   _SubDomainIt ite = _relevSubDomainsIntra[type].end();
   
   int i, tag=0, start = 0;
   for (i=0, it=itb; it!=ite; ++it, ++i) {
      int sbufsize = (it->second)->sendMsgSize();
      int number = it->first;
      assert(recvBuf != NULL);
      
      MPI_Irecv((void *) &(recvBuf[start]), sbufsize, (MPI_Datatype) MPI_BYTE, 
              number, tag, _pComm, &request[i]);
      
      start += sbufsize;
      assert (start<=myDomain->recvBufferSize());
   }
   
   // make sure that all the receives are posted before sending
   // think about defining a new communicator and synchronizing it only
   MPI_Barrier(_pComm);
   
   for (i=0, it=itb; it!=ite; ++it, ++i) {
      MPI_Rsend((void *) myDomain->sendBuffer(), myDomain->sendMsgSize(), 
              (MPI_Datatype) MPI_BYTE, it->first, tag, _pComm);
   }
   
   if (_relevSubDomainsIntra[type].size()) {
      MPI_Waitall(_relevSubDomainsIntra[type].size(), request, status);
   }
   
   if (request) std::free(request);
   if (status) std::free(status);
}

// helper - for debugging purposes
template <typename V>
template <typename Iterator>
void DomainCoupler<V>::relevantSubDomainNumbers(Iterator domainNumbers, int t, char *tch) {
   std::string commT(tch);
   _SubDomainIt it, ite;
   
   if (commT== "Inter") {
      it = _relevSubDomainsInter[t].begin();
      ite = _relevSubDomainsInter[t].end();  
   } else if (commT=="Intra") {
      it = _relevSubDomainsIntra[t].begin();
      ite = _relevSubDomainsIntra[t].end();
   } else 
      return;
   
   for ( ; it!=ite; ++it) *(domainNumbers++) = it->first;
   
   return;
}

template <typename V>
template <typename T>
void DomainCoupler<V>::registerSendBuffer(const SolvType& type, T* const start, int size) {
   int rank;
   MPI_Comm_rank(_pComm, &rank);
   
   _SubDomainIt it, ite=_subDomains[type].end();
   
   it = _subDomains[type].find(rank);
   if (it==ite) return;
   
   SubDomain<V>* myDomain = it->second;
   
   size_t unit = sizeof(T);
   size *= unit;
   
   char *startCh = reinterpret_cast<char *>(start);
   
   myDomain->registerSendBuffer(startCh, size);
   
   return;
}

template <typename V>
template <typename RT>
RT* DomainCoupler<V>::recvBuffer(const SolvType& type, const RT& notUsed) {
   // notUsed is only "used" by the compiler
   char *tmp=NULL;
   int rank;
   MPI_Comm_rank(_pComm, &rank);
   
   _SubDomainIt it, ite=_subDomains[type].end();
   it = _subDomains[type].find(rank);
   if (it==ite) return NULL;
   SubDomain<V>* myDomain = it->second;
   
   tmp = myDomain->recvBuffer();
   
   RT *tmpRT = reinterpret_cast<RT *>(tmp);
   
   return tmpRT;
}

template <typename V>
template <typename T>
int DomainCoupler<V>::recvBufferSize(const SolvType& type) {
   int tmp=0, rank;
   MPI_Comm_rank(_pComm, &rank);
   
   _SubDomainIt it, ite=_subDomains[type].end();
   it = _subDomains[type].find(rank);
   if (it==ite) return tmp;
   SubDomain<V>* myDomain = it->second;
   
   tmp = myDomain->recvBufferSize();
   
   size_t unit = sizeof(T);
   tmp /= unit;
   
   return tmp;
}
   
} // namespace pico

} // namespace DAPTA

---