src/ParaMEDMEM/ExplicitCoincidentDEC.cxx

   1 // Copyright (C) 2007-2020  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include <mpi.h>
  21 #include "CommInterface.hxx"
  22 #include "Topology.hxx"
  23 #include "BlockTopology.hxx"
  24 #include "ComponentTopology.hxx"
  25 #include "ParaFIELD.hxx"
  26 #include "MPIProcessorGroup.hxx"
  27 #include "ExplicitCoincidentDEC.hxx"
  28 #include "ExplicitMapping.hxx"
  29 #include "InterpKernelUtilities.hxx"
  30
  31 using namespace std;
  32
  33 namespace MEDCoupling
  34 {
  35
  36   /*!
  37     \anchor ExplicitCoincidentDEC-det
  38     \class ExplicitCoincidentDEC
  39
  40
  41     This class aims at \ref interpolation "remapping fields" that have identical
  42     supports (=the same underlying mesh) but different parallel topologies
  43     (=different sub-domains in the mesh). It can be used to couple
  44     together multi-physics codes that operate on the same domain
  45     with different partitioning.
  46
  47     It is very similar to what the \ref StructuredCoincidentDEC-det "StructuredCoincidentDEC"
  48     does, except that it works with an arbitrary user-defined topology.
  49
  50     The remapping between the two supports is based on identity of global
  51     ids, instead of geometrical considerations (as it is the case for
  52     \ref InterpKernelDEC-det "InterpKernelDEC").
  53     Therefore, beware that this \ref para-dec "DEC" can not be used
  54     for coincident meshes if they do *not* have the exact same numbering.
  55
  56     With this \ref para-dec "DEC" no projection, and no interpolation of the field data is done, contrary
  57     to what happens in \ref InterpKernelDEC-det "InterpKernelDEC". It is just
  58     a matter of allocating the values from one side to the other, using directly the cell
  59     identifiers.
  60
  61     As all the other DECs, its usage requires two phases :
  62     - a setup phase during which the topologies are exchanged so that
  63     the target side knows from which processors it should expect
  64     the data.
  65     - a send/recv phase during which the field data is actually transferred.
  66
  67     This example illustrates the sending of a field with
  68     the \c ExplicitCoincidentDEC :
  69     \code
  70     ...
  71     ExplicitCoincidentDEC dec(groupA, groupB);
  72     dec.attachLocalField(field);
  73     dec.synchronize();
  74     if (groupA.containsMyRank())
  75     dec.recvData();
  76     else if (groupB.containsMyRank())
  77     dec.sendData();
  78     ...
  79     \endcode
  80
  81     Creating a ParaFIELD to be attached to the %DEC is done in exactly the same way as for
  82     the other DECs, if only the partitioning of the support mesh differs.
  83     In the case where the
  84     fields have also different *component* topologies, creating the ParaFIELD
  85     requires some more effort. See the \ref para-over "parallelism" section for more details.
  86   */
  87
  88
  89   /*! Constructor
  90    */
  91   ExplicitCoincidentDEC::ExplicitCoincidentDEC():
  92       _toposource(0),_topotarget(0),
  93       _targetgroup(0), _sourcegroup(0),
  94       _sendcounts(0), _recvcounts(0),
  95       _senddispls(0), _recvdispls(0),
  96       _recvbuffer(0), _sendbuffer(0),
  97       _distant_elems(), _explicit_mapping()
  98   {
  99   }
 100
 101   ExplicitCoincidentDEC::~ExplicitCoincidentDEC()
 102   {
 103   }
 104
 105   /*! Synchronization process for exchanging topologies
 106    */
 107   void ExplicitCoincidentDEC::synchronize()
 108   {
 109     if (_source_group->containsMyRank())
 110       {
 111         _toposource = dynamic_cast<ExplicitTopology*>(_local_field->getTopology());
 112         _sourcegroup= _toposource->getProcGroup()->createProcGroup();
 113         _targetgroup=_toposource->getProcGroup()->createComplementProcGroup();
 114       }
 115     if (_target_group->containsMyRank())
 116       {
 117         _topotarget = dynamic_cast<ExplicitTopology*>(_local_field->getTopology());
 118         _sourcegroup= _topotarget->getProcGroup()->createComplementProcGroup();
 119         _targetgroup=_topotarget->getProcGroup()->createProcGroup();
 120       }
 121
 122     // Exchanging
 123
 124     // Transmitting source topology to target code
 125     broadcastTopology(_toposource,_topotarget,1000);
 126     transferMappingToSource();
 127   }
 128
 129   /*! Creates the arrays necessary for the data transfer
 130    * and fills the send array with the values of the
 131    * source field
 132    *  */
 133   void ExplicitCoincidentDEC::prepareSourceDE()
 134   {
 135     ////////////////////////////////////
 136     //Step 1 : buffer array creation
 137
 138     if (!_toposource->getProcGroup()->containsMyRank())
 139       return;
 140     MPIProcessorGroup* group=new MPIProcessorGroup(_sourcegroup->getCommInterface());
 141
 142     // Warning : the size of the target side is implicitly deduced
 143     //from the size of MPI_COMM_WORLD
 144     int target_size = _toposource->getProcGroup()->getCommInterface().worldSize()- _toposource->getProcGroup()->size()  ;
 145
 146     vector<int>* target_arrays=new vector<int>[target_size];
 147
 148     mcIdType nb_local = _toposource-> getNbLocalElements();
 149
 150     std::size_t union_size=group->size();
 151
 152     _sendcounts=new int[union_size];
 153     _senddispls=new int[union_size];
 154     _recvcounts=new int[union_size];
 155     _recvdispls=new int[union_size];
 156
 157     for (std::size_t i=0; i< union_size; i++)
 158       {
 159         _sendcounts[i]=0;
 160         _recvcounts[i]=0;
 161         _recvdispls[i]=0;
 162       }
 163     _senddispls[0]=0;
 164
 165     int* counts=_explicit_mapping.getCounts();
 166     for (int i=0; i<group->size(); i++)
 167       _sendcounts[i]=counts[i];
 168
 169     for (int iproc=1; iproc<group->size();iproc++)
 170       _senddispls[iproc]=_senddispls[iproc-1]+_sendcounts[iproc-1];
 171
 172     _sendbuffer = new double [nb_local * _toposource->getNbComponents()];
 173
 174     /////////////////////////////////////////////////////////////
 175     //Step 2 : filling the buffers with the source field values
 176
 177     int* counter=new int [target_size];
 178     counter[0]=0;
 179     for (int i=1; i<target_size; i++)
 180       counter[i]=counter[i-1]+(int)target_arrays[i-1].size();
 181
 182
 183     const double* value = _local_field->getField()->getArray()->getPointer();
 184
 185     int* bufferindex= _explicit_mapping.getBufferIndex();
 186
 187     for (int ielem=0; ielem<nb_local; ielem++)
 188       {
 189         int ncomp = _toposource->getNbComponents();
 190         for (int icomp=0; icomp<ncomp; icomp++)
 191           {
 192             _sendbuffer[ielem*ncomp+icomp]=value[bufferindex[ielem]*ncomp+icomp];
 193           }
 194       }
 195     delete[] target_arrays;
 196     delete[] counter;
 197   }
 198
 199   /*!
 200    *  Creates the buffers for receiving the fields on the target side
 201    */
 202   void ExplicitCoincidentDEC::prepareTargetDE()
 203   {
 204     if (!_topotarget->getProcGroup()->containsMyRank())
 205       return;
 206     MPIProcessorGroup* group=new MPIProcessorGroup(_topotarget->getProcGroup()->getCommInterface());
 207
 208     vector < vector <mcIdType> > source_arrays(_sourcegroup->size());
 209     mcIdType nb_local = _topotarget-> getNbLocalElements();
 210     for (mcIdType ielem=0; ielem< nb_local ; ielem++)
 211       {
 212         //pair<int,mcIdType> source_local =_distant_elems[ielem];
 213         pair <int,mcIdType> source_local=_explicit_mapping.getDistantNumbering(ielem);
 214         source_arrays[source_local.first].push_back(source_local.second);
 215       }
 216     std::size_t union_size=group->size();
 217     _recvcounts=new int[union_size];
 218     _recvdispls=new int[union_size];
 219     _sendcounts=new int[union_size];
 220     _senddispls=new int[union_size];
 221
 222     for (std::size_t i=0; i< union_size; i++)
 223       {
 224         _sendcounts[i]=0;
 225         _recvcounts[i]=0;
 226         _recvdispls[i]=0;
 227       }
 228     for (int iproc=0; iproc < _sourcegroup->size(); iproc++)
 229       {
 230         //converts the rank in target to the rank in union communicator
 231         int unionrank=group->translateRank(_sourcegroup,iproc);
 232         _recvcounts[unionrank]=(int)(source_arrays[iproc].size()*_topotarget->getNbComponents());
 233       }
 234     for (std::size_t i=1; i<union_size; i++)
 235       _recvdispls[i]=_recvdispls[i-1]+_recvcounts[i-1];
 236     _recvbuffer=new double[nb_local*_topotarget->getNbComponents()];
 237
 238   }
 239
 240
 241   /*!
 242    * Synchronizing a topology so that all the groups get it.
 243    *
 244    * \param toposend Topology that is transmitted. It is read on processes where it already exists, and it is created and filled on others.
 245    * \param toporecv Topology which is received.
 246    * \param tag Communication tag associated with this operation.
 247    */
 248   void ExplicitCoincidentDEC::broadcastTopology(const ExplicitTopology* toposend, ExplicitTopology* toporecv, int tag)
 249   {
 250     MPI_Status status;
 251
 252     mcIdType* serializer=0;
 253     mcIdType size;
 254
 255     MPIProcessorGroup* group=new MPIProcessorGroup(*_comm_interface);
 256
 257     // The send processors serialize the send topology
 258     // and send the buffers to the recv procs
 259     if (toposend !=0 && toposend->getProcGroup()->containsMyRank())
 260       {
 261         toposend->serialize(serializer, size);
 262         for (int iproc=0; iproc< group->size(); iproc++)
 263           {
 264             int itarget=iproc;
 265             if (!toposend->getProcGroup()->contains(itarget))
 266               {
 267                 _comm_interface->send(&size,1,MPI_ID_TYPE, itarget,tag+itarget,*(group->getComm()));
 268                 _comm_interface->send(serializer, (int)size, MPI_ID_TYPE, itarget, tag+itarget,*(group->getComm()));
 269               }
 270           }
 271       }
 272     else
 273       {
 274         vector <int> size2(group->size());
 275         int myworldrank=group->myRank();
 276         for (int iproc=0; iproc<group->size();iproc++)
 277           {
 278             int isource = iproc;
 279             if (!toporecv->getProcGroup()->contains(isource))
 280               {
 281                 mcIdType nbelem;
 282                 _comm_interface->recv(&nbelem, 1, MPI_ID_TYPE, isource, tag+myworldrank, *(group->getComm()), &status);
 283                 mcIdType* buffer = new mcIdType[nbelem];
 284                 _comm_interface->recv(buffer, (int)nbelem, MPI_ID_TYPE, isource,tag+myworldrank, *(group->getComm()), &status);
 285
 286                 ExplicitTopology* topotemp=new ExplicitTopology();
 287                 topotemp->unserialize(buffer, *_comm_interface);
 288                 delete[] buffer;
 289
 290                 for (mcIdType ielem=0; ielem<toporecv->getNbLocalElements(); ielem++)
 291                   {
 292                     mcIdType global = toporecv->localToGlobal(ielem);
 293                     mcIdType sendlocal=topotemp->globalToLocal(global);
 294                     if (sendlocal!=-1)
 295                       {
 296                         size2[iproc]++;
 297                         _explicit_mapping.pushBackElem(make_pair(iproc,sendlocal));
 298                       }
 299                   }
 300                 delete topotemp;
 301               }
 302           }
 303       }
 304     MESSAGE (" rank "<<group->myRank()<< " broadcastTopology is over");
 305   }
 306
 307   void ExplicitCoincidentDEC::transferMappingToSource()
 308   {
 309
 310     MPIProcessorGroup* group=new MPIProcessorGroup(*_comm_interface);
 311
 312     // sending source->target mapping which is stored by target
 313     //in _distant_elems from target to source
 314     if (_topotarget!=0 && _topotarget->getProcGroup()->containsMyRank())
 315       {
 316         int world_size = _topotarget->getProcGroup()->getCommInterface().worldSize()  ;
 317         int* nb_transfer_union=new int[world_size];
 318         int* dummy_recv=new int[world_size];
 319         for (int i=0; i<world_size; i++)
 320           nb_transfer_union[i]=0;
 321         //converts the rank in target to the rank in union communicator
 322
 323         for (int i=0; i<  _explicit_mapping.nbDistantDomains(); i++)
 324           {
 325             int unionrank=group->translateRank(_sourcegroup,_explicit_mapping.getDistantDomain(i));
 326             nb_transfer_union[unionrank]=_explicit_mapping.getNbDistantElems(i);
 327           }
 328         _comm_interface->allToAll(nb_transfer_union, 1, MPI_INT, dummy_recv, 1, MPI_INT, MPI_COMM_WORLD);
 329
 330         mcIdType* sendbuffer= _explicit_mapping.serialize(_topotarget->getProcGroup()->myRank());
 331
 332         int* sendcounts= new int [world_size];
 333         int* senddispls = new int [world_size];
 334         for (int i=0; i< world_size; i++)
 335           {
 336             sendcounts[i]=2*nb_transfer_union[i];
 337             if (i==0)
 338               senddispls[i]=0;
 339             else
 340               senddispls[i]=senddispls[i-1]+sendcounts[i-1];
 341           }
 342         int* recvcounts=new int[world_size];
 343         int* recvdispls=new int[world_size];
 344         int *dummyrecv=0;
 345         for (int i=0; i <world_size; i++)
 346           {
 347             recvcounts[i]=0;
 348             recvdispls[i]=0;
 349           }
 350         _comm_interface->allToAllV(sendbuffer, sendcounts, senddispls, MPI_ID_TYPE, dummyrecv, recvcounts, senddispls, MPI_ID_TYPE, MPI_COMM_WORLD);
 351
 352       }
 353     //receiving in the source subdomains the mapping sent by targets
 354     else
 355       {
 356         int world_size = _toposource->getProcGroup()->getCommInterface().worldSize()  ;
 357         int* nb_transfer_union=new int[world_size];
 358         int* dummy_send=new int[world_size];
 359         for (int i=0; i<world_size; i++)
 360           dummy_send[i]=0;
 361         _comm_interface->allToAll(dummy_send, 1, MPI_INT, nb_transfer_union, 1, MPI_INT, MPI_COMM_WORLD);
 362
 363         int total_size=0;
 364         for (int i=0; i< world_size; i++)
 365           total_size+=nb_transfer_union[i];
 366         int nbtarget = _targetgroup->size();
 367         int* targetranks = new int[ nbtarget];
 368         for (int i=0; i<nbtarget; i++)
 369           targetranks[i]=group->translateRank(_targetgroup,i);
 370         mcIdType* mappingbuffer= new mcIdType [total_size*2];
 371         int* sendcounts= new int [world_size];
 372         int* senddispls = new int [world_size];
 373         int* recvcounts=new int[world_size];
 374         int* recvdispls=new int[world_size];
 375         for (int i=0; i< world_size; i++)
 376           {
 377             recvcounts[i]=2*nb_transfer_union[i];
 378             if (i==0)
 379               recvdispls[i]=0;
 380             else
 381               recvdispls[i]=recvdispls[i-1]+recvcounts[i-1];
 382           }
 383
 384         int *dummysend=0;
 385         for (int i=0; i <world_size; i++)
 386           {
 387             sendcounts[i]=0;
 388             senddispls[i]=0;
 389           }
 390         _comm_interface->allToAllV(dummysend, sendcounts, senddispls, MPI_ID_TYPE, mappingbuffer, recvcounts, recvdispls, MPI_ID_TYPE, MPI_COMM_WORLD);
 391         _explicit_mapping.unserialize(world_size,nb_transfer_union,nbtarget, targetranks, mappingbuffer);
 392       }
 393   }
 394
 395   void ExplicitCoincidentDEC::recvData()
 396   {
 397     //MPI_COMM_WORLD is used instead of group because there is no
 398     //mechanism for creating the union group yet
 399     MESSAGE("recvData");
 400
 401     cout<<"start AllToAll"<<endl;
 402     _comm_interface->allToAllV(_sendbuffer, _sendcounts, _senddispls, MPI_DOUBLE,
 403                                _recvbuffer, _recvcounts, _recvdispls, MPI_DOUBLE,MPI_COMM_WORLD);
 404     cout<<"end AllToAll"<<endl;
 405     mcIdType nb_local = _topotarget->getNbLocalElements();
 406     double* value=new double[nb_local*_topotarget->getNbComponents()];
 407
 408     vector<int> counters(_sourcegroup->size());
 409     counters[0]=0;
 410     for (int i=0; i<_sourcegroup->size()-1; i++)
 411       {
 412         MPIProcessorGroup* group=new MPIProcessorGroup(*_comm_interface);
 413         int worldrank=group->translateRank(_sourcegroup,i);
 414         counters[i+1]=counters[i]+_recvcounts[worldrank];
 415       }
 416
 417     for (int ielem=0; ielem<nb_local ; ielem++)
 418       {
 419         pair<int,int> distant_numbering=_explicit_mapping.getDistantNumbering(ielem);
 420         int iproc=distant_numbering.first;
 421         int ncomp =  _topotarget->getNbComponents();
 422         for (int icomp=0; icomp< ncomp; icomp++)
 423           value[ielem*ncomp+icomp]=_recvbuffer[counters[iproc]*ncomp+icomp];
 424         counters[iproc]++;
 425       }
 426     _local_field->getField()->getArray()->useArray(value,true,DeallocType::CPP_DEALLOC,nb_local,_topotarget->getNbComponents());
 427   }
 428
 429   void ExplicitCoincidentDEC::sendData()
 430   {
 431     MESSAGE ("sendData");
 432     for (int i=0; i< 4; i++)
 433       cout << _sendcounts[i]<<" ";
 434     cout <<endl;
 435     for (int i=0; i< 4; i++)
 436       cout << _senddispls[i]<<" ";
 437     cout <<endl;
 438     //MPI_COMM_WORLD is used instead of group because there is no
 439     //mechanism for creating the union group yet
 440     cout <<"start AllToAll"<<endl;
 441     _comm_interface->allToAllV(_sendbuffer, _sendcounts, _senddispls, MPI_DOUBLE,
 442                                _recvbuffer, _recvcounts, _recvdispls, MPI_DOUBLE,MPI_COMM_WORLD);
 443   }
 444 }
 445