src/ParaMEDMEM/StructuredCoincidentDEC.cxx

   1 //  Copyright (C) 2007-2008  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.
   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 #include <mpi.h>
  20 #include "CommInterface.hxx"
  21 #include "Topology.hxx"
  22 #include "BlockTopology.hxx"
  23 #include "ComponentTopology.hxx"
  24 #include "ParaFIELD.hxx"
  25 #include "MPIProcessorGroup.hxx"
  26 #include "StructuredCoincidentDEC.hxx"
  27 #include "InterpKernelUtilities.hxx"
  28
  29 #include <iostream>
  30
  31 using namespace std;
  32
  33 namespace ParaMEDMEM
  34 {
  35
  36   /*! \defgroup structuredcoincidentdec StructuredCoincidentDEC
  37
  38     This class is meant for remapping fields that have identical
  39     supports with different parallel topologies. It can be used to couple
  40     together multiphysics codes that operate on the same domain
  41     with different partitionings, which can be useful if one of
  42     the computation is much faster than the other. It can also be used
  43     to couple together codes that share an interface that was generated
  44     in the same manner (with identical global ids).
  45     Also, this DEC can be used for fields that have component topologies,
  46     i.e., components that are scattered over several processors.
  47
  48     The remapping between the two supports is based on identity of global
  49     ids, instead of geometrical considerations as it is the case for
  50     NonCoincidentDEC and IntersectionDEC. Therefore, this DEC must not be used
  51     for coincident meshes that do not have the same numbering.
  52
  53     As all the other DECs, its use is made of two phases :
  54     - a setup phase during whih the topologies are exchanged so that
  55     the target side knows from which processors it should expect
  56     the data.
  57     - a send/recv phase during which the field data is actually transferred.
  58
  59     This example illustrates the sending of a field with
  60     the DEC :
  61     \code
  62     ...
  63     StructuredCoincidentDEC dec(groupA, groupB);
  64     dec.attachLocalField(field);
  65     dec.synchronize();
  66     if (groupA.containsMyRank())
  67     dec.recvData();
  68     else if (groupB.containsMyRank())
  69     dec.sendData();
  70     ...
  71     \endcode
  72
  73     Creating a ParaFIELD to be attached to the DEC is exactly the same as for
  74     other DECs in the case when the remapping concerns similar meshes
  75     that only have different partitionings. In the case when the
  76     fields have also different component topologies, creating the ParaFIELD
  77     requires some more effort. See \ref parafield section for more details.
  78   */
  79
  80
  81   StructuredCoincidentDEC::StructuredCoincidentDEC():_topo_source(0),_topo_target(0),
  82                                                      _recv_buffer(0),_send_buffer(0),
  83                                                      _recv_counts(0),_send_counts(0),
  84                                                      _recv_displs(0),_send_displs(0)
  85   {
  86   }
  87
  88
  89   StructuredCoincidentDEC::~StructuredCoincidentDEC()
  90   {
  91     delete[] _send_buffer;
  92     delete[] _recv_buffer;
  93     delete[]_send_displs;
  94     delete[] _recv_displs;
  95     delete[] _send_counts;
  96     delete[] _recv_counts;
  97     if (! _source_group->containsMyRank())
  98       delete _topo_source;
  99     if(!_target_group->containsMyRank())
 100       delete _topo_target;
 101   }
 102
 103   /*!
 104     \addtogroup structuredcoincidentdec
 105     @{
 106   */
 107   StructuredCoincidentDEC::StructuredCoincidentDEC(ProcessorGroup& local_group, ProcessorGroup& distant_group):DEC(local_group,distant_group),_topo_source(0),_topo_target(0),_recv_buffer(0),_send_buffer(0)
 108   {
 109   }
 110
 111   /*! Synchronization process for exchanging topologies
 112    */
 113   void StructuredCoincidentDEC::synchronizeTopology()
 114   {
 115     if (_source_group->containsMyRank())
 116       _topo_source = dynamic_cast<BlockTopology*>(_local_field->getTopology());
 117     if (_target_group->containsMyRank())
 118       _topo_target = dynamic_cast<BlockTopology*>(_local_field->getTopology());
 119
 120     // Transmitting source topology to target code
 121     broadcastTopology(_topo_source,1000);
 122     // Transmitting target topology to source code
 123     broadcastTopology(_topo_target,2000);
 124     if (_topo_source->getNbElements() != _topo_target->getNbElements())
 125       throw INTERP_KERNEL::Exception("Incompatible dimensions for target and source topologies");
 126
 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 StructuredCoincidentDEC::prepareSourceDE()
 134   {
 135     ////////////////////////////////////
 136     //Step 1 : _buffer array creation
 137
 138     if (!_topo_source->getProcGroup()->containsMyRank())
 139       return;
 140     MPIProcessorGroup* group=new MPIProcessorGroup(_topo_source->getProcGroup()->getCommInterface());
 141
 142     int myranksource = _topo_source->getProcGroup()->myRank();
 143
 144     vector <int>* target_arrays=new vector<int>[_topo_target->getProcGroup()->size()];
 145
 146     //cout<<" topotarget size"<<  _topo_target->getProcGroup()->size()<<endl;
 147
 148     int nb_local = _topo_source-> getNbLocalElements();
 149     for (int ielem=0; ielem< nb_local ; ielem++)
 150       {
 151         //  cout <<"source local :"<<myranksource<<","<<ielem<<endl;
 152         int global = _topo_source->localToGlobal(make_pair(myranksource, ielem));
 153         //  cout << "global "<<global<<endl;
 154         pair<int,int> target_local =_topo_target->globalToLocal(global);
 155         //  cout << "target local : "<<target_local.first<<","<<target_local.second<<endl;
 156         target_arrays[target_local.first].push_back(target_local.second);
 157       }
 158
 159     int union_size=group->size();
 160
 161     _send_counts=new int[union_size];
 162     _send_displs=new int[union_size];
 163     _recv_counts=new int[union_size];
 164     _recv_displs=new int[union_size];
 165
 166     for (int i=0; i< union_size; i++)
 167       {
 168         _send_counts[i]=0;
 169         _recv_counts[i]=0;
 170         _recv_displs[i]=0;
 171       }
 172     _send_displs[0]=0;
 173
 174     for (int iproc=0; iproc < _topo_target->getProcGroup()->size(); iproc++)
 175       {
 176         //converts the rank in target to the rank in union communicator
 177         int unionrank=group->translateRank(_topo_target->getProcGroup(),iproc);
 178         _send_counts[unionrank]=target_arrays[iproc].size();
 179       }
 180
 181     for (int iproc=1; iproc<group->size();iproc++)
 182       _send_displs[iproc]=_send_displs[iproc-1]+_send_counts[iproc-1];
 183
 184     _send_buffer = new double [nb_local ];
 185
 186     /////////////////////////////////////////////////////////////
 187     //Step 2 : filling the _buffers with the source field values
 188
 189     int* counter=new int [_topo_target->getProcGroup()->size()];
 190     counter[0]=0;
 191     for (int i=1; i<_topo_target->getProcGroup()->size(); i++)
 192       counter[i]=counter[i-1]+target_arrays[i-1].size();
 193
 194
 195     const double* value = _local_field->getField()->getArray()->getPointer();
 196     //cout << "Nb local " << nb_local<<endl;
 197     for (int ielem=0; ielem<nb_local ; ielem++)
 198       {
 199         int global = _topo_source->localToGlobal(make_pair(myranksource, ielem));
 200         pair<int,int> target_local =_topo_target->globalToLocal(global);
 201         //cout <<"global : "<< global<<" local :"<<target_local.first<<" "<<target_local.second;
 202         //cout <<"counter[]"<<counter[target_local.first]<<endl;
 203         _send_buffer[counter[target_local.first]++]=value[ielem];
 204
 205       }
 206     delete[] target_arrays;
 207     delete[] counter;
 208     delete group;
 209   }
 210
 211   /*!
 212    *  Creates the _buffers for receiving the fields on the target side
 213    */
 214   void StructuredCoincidentDEC::prepareTargetDE()
 215   {
 216     if (!_topo_target->getProcGroup()->containsMyRank())
 217       return;
 218     MPIProcessorGroup* group=new MPIProcessorGroup(_topo_source->getProcGroup()->getCommInterface());
 219
 220     int myranktarget = _topo_target->getProcGroup()->myRank();
 221
 222     vector < vector <int> > source_arrays(_topo_source->getProcGroup()->size());
 223     int nb_local = _topo_target-> getNbLocalElements();
 224     for (int ielem=0; ielem< nb_local ; ielem++)
 225       {
 226         //  cout <<"TS target local :"<<myranktarget<<","<<ielem<<endl;
 227         int global = _topo_target->localToGlobal(make_pair(myranktarget, ielem));
 228         //cout << "TS global "<<global<<endl;
 229         pair<int,int> source_local =_topo_source->globalToLocal(global);
 230         //  cout << "TS source local : "<<source_local.first<<","<<source_local.second<<endl;
 231         source_arrays[source_local.first].push_back(source_local.second);
 232       }
 233     int union_size=group->size();
 234     _recv_counts=new int[union_size];
 235     _recv_displs=new int[union_size];
 236     _send_counts=new int[union_size];
 237     _send_displs=new int[union_size];
 238
 239     for (int i=0; i< union_size; i++)
 240       {
 241         _send_counts[i]=0;
 242         _recv_counts[i]=0;
 243         _recv_displs[i]=0;
 244       }
 245     for (int iproc=0; iproc < _topo_source->getProcGroup()->size(); iproc++)
 246       {
 247         //converts the rank in target to the rank in union communicator
 248         int unionrank=group->translateRank(_topo_source->getProcGroup(),iproc);
 249         _recv_counts[unionrank]=source_arrays[iproc].size();
 250       }
 251     for (int i=1; i<union_size; i++)
 252       _recv_displs[i]=_recv_displs[i-1]+_recv_counts[i-1];
 253     _recv_buffer=new double[nb_local];
 254
 255     delete group;
 256   }
 257
 258
 259   /*!
 260    * Synchronizing a topology so that all the
 261    * group possesses it.
 262    *
 263    * \param topo Topology that is transmitted. It is read on processes where it already exists, and it is created and filled on others.
 264    * \param tag Communication tag associated with this operation.
 265    */
 266   void StructuredCoincidentDEC::broadcastTopology(BlockTopology*& topo, int tag)
 267   {
 268     MPI_Status status;
 269
 270     int* serializer=0;
 271     int size;
 272
 273     MPIProcessorGroup* group=new MPIProcessorGroup(*_comm_interface);
 274
 275     // The master proc creates a send buffer containing
 276     // a serialized topology
 277     int rank_master;
 278
 279     if (topo!=0 && topo->getProcGroup()->myRank()==0)
 280       {
 281         MESSAGE ("Master rank");
 282         topo->serialize(serializer, size);
 283         rank_master = group->translateRank(topo->getProcGroup(),0);
 284         MESSAGE("Master rank world number is "<<rank_master);
 285         MESSAGE("World Size is "<<group->size());
 286         for (int i=0; i< group->size(); i++)
 287           {
 288             if (i!= rank_master)
 289               _comm_interface->send(&rank_master,1,MPI_INT, i,tag+i,*(group->getComm()));
 290           }
 291       }
 292     else
 293       {
 294         MESSAGE(" rank "<<group->myRank()<< " waiting ...");
 295         _comm_interface->recv(&rank_master, 1,MPI_INT, MPI_ANY_SOURCE, tag+group->myRank(), *(group->getComm()),&status);
 296         MESSAGE(" rank "<<group->myRank()<< "received master rank"<<rank_master);
 297       }
 298     // The topology is broadcasted to all processsors in the group
 299     _comm_interface->broadcast(&size, 1,MPI_INT,rank_master,*(group->getComm()));
 300
 301     int* buffer=new int[size];
 302     if (topo!=0 && topo->getProcGroup()->myRank()==0)
 303       copy(serializer, serializer+size, buffer);
 304     _comm_interface->broadcast(buffer,size,MPI_INT,rank_master,*(group->getComm()));
 305
 306     // Processors which did not possess the source topology
 307     // unserialize it
 308
 309     BlockTopology* topotemp=new BlockTopology();
 310     topotemp->unserialize(buffer, *_comm_interface);
 311
 312     if (topo==0)
 313       topo=topotemp;
 314     else
 315       delete topotemp;
 316
 317     // Memory cleaning
 318     delete[] buffer;
 319     if (serializer!=0)
 320       delete[] serializer;
 321     MESSAGE (" rank "<<group->myRank()<< " unserialize is over");
 322     delete group;
 323   }
 324
 325
 326
 327   void StructuredCoincidentDEC::recvData()
 328   {
 329     //MPI_COMM_WORLD is used instead of group because there is no
 330     //mechanism for creating the union group yet
 331     MESSAGE("recvData");
 332     for (int i=0; i< 4; i++)
 333       cout << _recv_counts[i]<<" ";
 334     cout <<endl;
 335     for (int i=0; i< 4; i++)
 336       cout << _recv_displs[i]<<" ";
 337     cout <<endl;
 338
 339     cout<<"start AllToAll"<<endl;
 340     MPI_Comm comm = *(dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
 341     _comm_interface->allToAllV(_send_buffer, _send_counts, _send_displs, MPI_DOUBLE,
 342                                _recv_buffer, _recv_counts, _recv_displs, MPI_DOUBLE,comm);
 343     cout<<"end AllToAll"<<endl;
 344
 345     int nb_local = _topo_target->getNbLocalElements();
 346     //double* value=new double[nb_local];
 347     double* value=const_cast<double*>(_local_field->getField()->getArray()->getPointer());
 348
 349     int myranktarget=_topo_target->getProcGroup()->myRank();
 350     vector<int> counters(_topo_source->getProcGroup()->size());
 351     counters[0]=0;
 352     for (int i=0; i<_topo_source->getProcGroup()->size()-1; i++)
 353       {
 354         MPIProcessorGroup* group=new MPIProcessorGroup(*_comm_interface);
 355         int worldrank=group->translateRank(_topo_source->getProcGroup(),i);
 356         counters[i+1]=counters[i]+_recv_counts[worldrank];
 357         delete group;
 358       }
 359
 360     for (int ielem=0; ielem<nb_local ; ielem++)
 361       {
 362         int global = _topo_target->localToGlobal(make_pair(myranktarget, ielem));
 363         pair<int,int> source_local =_topo_source->globalToLocal(global);
 364         value[ielem]=_recv_buffer[counters[source_local.first]++];
 365       }
 366
 367
 368     //_local_field->getField()->setValue(value);
 369   }
 370
 371   void StructuredCoincidentDEC::sendData()
 372   {
 373     MESSAGE ("sendData");
 374     for (int i=0; i< 4; i++)
 375       cout << _send_counts[i]<<" ";
 376     cout <<endl;
 377     for (int i=0; i< 4; i++)
 378       cout << _send_displs[i]<<" ";
 379     cout <<endl;
 380     cout <<"start AllToAll"<<endl;
 381     MPI_Comm comm = *(dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
 382     _comm_interface->allToAllV(_send_buffer, _send_counts, _send_displs, MPI_DOUBLE,
 383                                _recv_buffer, _recv_counts, _recv_displs, MPI_DOUBLE,comm);
 384     cout<<"end AllToAll"<<endl;
 385   }
 386
 387   /*! Prepares a DEC for data exchange
 388
 389     This method broadcasts the topologies from source to target
 390     so that the target side can analyse from which processors it
 391     is expected to receive data.
 392   */
 393
 394   void StructuredCoincidentDEC::synchronize()
 395   {
 396     if (_source_group->containsMyRank())
 397       {
 398         synchronizeTopology();
 399         prepareSourceDE();
 400       }
 401     else if (_target_group->containsMyRank())
 402       {
 403         synchronizeTopology();
 404         prepareTargetDE();
 405       }
 406   }
 407   /*!
 408     @}
 409   */
 410 }
 411