src/ParaMEDMEM/BlockTopology.cxx

   1 // Copyright (C) 2007-2015  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 "BlockTopology.hxx"
  21 #include "MEDCouplingMemArray.hxx"
  22 #include "MEDCouplingCMesh.hxx"
  23 #include "CommInterface.hxx"
  24 #include "ProcessorGroup.hxx"
  25 #include "MPIProcessorGroup.hxx"
  26 #include "ComponentTopology.hxx"
  27 #include "InterpKernelUtilities.hxx"
  28
  29 #include <vector>
  30 #include <algorithm>
  31 #include <utility>
  32 #include <iostream>
  33
  34 using namespace std;
  35
  36 namespace ParaMEDMEM
  37 {
  38
  39   //!converts a pair <subdomainid,local> to a global number
  40   std::pair<int,int> BlockTopology::globalToLocal(const int global) const
  41   {
  42     int subdomain_id=0;
  43     int position=global;
  44     int size=_nb_elems;
  45     int size_procs=_proc_group->size();
  46     int increment=size;
  47     vector<int>axis_position(_dimension);
  48     vector<int>axis_offset(_dimension);
  49     for (int idim=0; idim<_dimension; idim++)
  50       {
  51         int axis_size=_local_array_indices[idim].size()-1;
  52         int axis_nb_elem=_local_array_indices[idim][axis_size];
  53         increment=increment/axis_nb_elem;
  54         int proc_increment = size_procs/(axis_size);
  55         int axis_pos=position/increment;
  56         position=position%increment;
  57         int iaxis=1;
  58         while (_local_array_indices[idim][iaxis]<=axis_pos)
  59           {
  60             subdomain_id+=proc_increment;
  61             iaxis++;
  62           }
  63         axis_position[idim]=axis_pos-_local_array_indices[idim][iaxis-1];
  64         axis_offset[idim]=iaxis;
  65       }
  66     int local=0;
  67     int local_increment=1;
  68     for (int idim=_dimension-1; idim>=0; idim--)
  69       {
  70         local+=axis_position[idim]*local_increment;
  71         local_increment*=_local_array_indices[idim][axis_offset[idim]]-_local_array_indices[idim][axis_offset[idim]-1];
  72       }
  73     return make_pair(subdomain_id,local);
  74   }
  75
  76   //!converts local number to a global number
  77   int BlockTopology::localToGlobal(const pair<int,int> local) const
  78   {
  79
  80     int subdomain_id=local.first;
  81     int global=0;
  82     int loc=local.second;
  83     int increment=_nb_elems;
  84     int proc_increment=_proc_group->size();
  85     int local_increment=getNbLocalElements();
  86     for (int idim=0; idim < _dimension; idim++)
  87       {
  88         int axis_size=_local_array_indices[idim].size()-1;
  89         int axis_nb_elem=_local_array_indices[idim][axis_size];
  90         increment=axis_nb_elem==0?0:increment/axis_nb_elem;
  91         proc_increment = proc_increment/(axis_size);
  92         int proc_axis=subdomain_id/proc_increment;
  93         subdomain_id=subdomain_id%proc_increment;
  94         int local_axis_nb_elem=_local_array_indices[idim][proc_axis+1]-_local_array_indices[idim][proc_axis];
  95         local_increment = (local_axis_nb_elem==0)?0:(local_increment/local_axis_nb_elem);
  96         int iaxis=((local_increment==0)?0:(loc/local_increment))+_local_array_indices[idim][proc_axis];
  97         global+=increment*iaxis;
  98         loc = (local_increment==0)?0:(loc%local_increment);
  99       }
 100     return global;
 101   }
 102
 103   //Retrieves the local number of elements
 104   int BlockTopology::getNbLocalElements()const
 105   {
 106     int position=_proc_group->myRank();
 107     int nb_elem = 1;
 108     int increment=1;
 109     for (int i=_dimension-1; i>=0; i--)
 110       {
 111         increment *=_nb_procs_per_dim[i];
 112         int idim=position%increment;
 113         position=position/increment;
 114         int imin=_local_array_indices[i][idim];
 115         int imax=_local_array_indices[i][idim+1];
 116         nb_elem*=(imax-imin);
 117       }
 118     return nb_elem;
 119   }
 120
 121   /*!
 122    * Constructor of a block topology from a grid.
 123    * This preliminary version simply splits along the first axis
 124    * instead of making the best choice with respect to the
 125    * values of the different axes.
 126    */
 127   BlockTopology::BlockTopology(const ProcessorGroup& group, MEDCouplingCMesh *grid):
 128     _dimension(grid->getSpaceDimension()), _proc_group(&group), _owns_processor_group(false)
 129   {
 130     vector <int> axis_length(_dimension);
 131     _nb_elems=1;
 132     for (int idim=0; idim <_dimension; idim++)
 133       {
 134         DataArrayDouble *arr=grid->getCoordsAt(idim);
 135         axis_length[idim]=arr->getNbOfElems();
 136         _nb_elems*=axis_length[idim];
 137       }
 138     //default splitting along 1st dimension
 139     _local_array_indices.resize(_dimension);
 140     _nb_procs_per_dim.resize(_dimension);
 141
 142     _local_array_indices[0].resize(_proc_group->size()+1);
 143     _local_array_indices[0][0]=0;
 144     _nb_procs_per_dim[0]=_proc_group->size();
 145
 146     for (int i=1; i<=_proc_group->size(); i++)
 147       {
 148         _local_array_indices[0][i]=_local_array_indices[0][i-1]+
 149           axis_length[0]/_proc_group->size();
 150         if (i<= axis_length[0]%_proc_group->size())
 151           _local_array_indices[0][i]+=1;
 152       }
 153     for (int i=1; i<_dimension; i++)
 154       {
 155         _local_array_indices[i].resize(2);
 156         _local_array_indices[i][0]=0;
 157         _local_array_indices[i][1]=axis_length[i];
 158         _nb_procs_per_dim[i]=1;
 159       }
 160     _cycle_type.resize(_dimension);
 161     for (int i=0; i<_dimension; i++)
 162       _cycle_type[i]=ParaMEDMEM::Block;
 163   }
 164
 165   /*!
 166    * Creation of a block topology by composing
 167    * a geometrical topology and a component topology.
 168    * This constructor is intended for creating fields
 169    * for which the parallel distribution is made on the
 170    * components of the field rather than on the geometrical
 171    * partitioning of the underlying mesh.
 172    *
 173    */
 174   BlockTopology::BlockTopology(const BlockTopology& geom_topo, const ComponentTopology& comp_topo):_owns_processor_group(false)
 175   {
 176     // so far, the block topology can only be created if the proc group
 177     // is either on geom_topo or on comp_topo
 178     if (geom_topo.getProcGroup()->size()>1 && comp_topo.nbBlocks()>1)
 179       throw INTERP_KERNEL::Exception(LOCALIZED("BlockTopology cannot yet be constructed with both complex geo and components topology"));
 180
 181     if (comp_topo.nbComponents()==1)
 182       {
 183         *this=geom_topo;
 184         return;
 185       }
 186     else
 187       {
 188         _dimension = geom_topo.getDimension()+1;
 189         if (comp_topo.nbBlocks()>1)
 190           _proc_group=comp_topo.getProcGroup();
 191         else
 192           _proc_group=geom_topo.getProcGroup();
 193         _local_array_indices=geom_topo._local_array_indices;
 194         vector<int> comp_indices = *(comp_topo.getBlockIndices());
 195         _local_array_indices.push_back(comp_indices);
 196         _nb_procs_per_dim=geom_topo._nb_procs_per_dim;
 197         _nb_procs_per_dim.push_back(comp_topo.nbBlocks());
 198         _cycle_type=geom_topo._cycle_type;
 199         _cycle_type.push_back(Block);
 200         _nb_elems=geom_topo.getNbElements()*comp_topo.nbComponents();
 201       }
 202   }
 203
 204   /*! Constructor for creating a one-dimensional
 205    * topology from a processor group and a local
 206    * number of elements on each processor
 207    *
 208    * The function must be called only by the processors belonging
 209    * to group \a group. Calling it from a processor not belonging
 210    * to \a group will cause an MPI error, while calling from a subset
 211    * of \a group will result in a deadlock.
 212    */
 213   BlockTopology::BlockTopology(const ProcessorGroup& group, int nb_elem):_dimension(1),_proc_group(&group),_owns_processor_group(false)
 214   {
 215     int* nbelems_per_proc = new int[group.size()];
 216     const MPIProcessorGroup* mpi_group=dynamic_cast<const MPIProcessorGroup*>(_proc_group);
 217     const MPI_Comm* comm=mpi_group->getComm();
 218     int nbtemp=nb_elem;
 219     mpi_group->getCommInterface().allGather(&nbtemp, 1, MPI_INT,
 220                                             nbelems_per_proc, 1, MPI_INT,
 221                                             *comm);
 222     _nb_elems=0;
 223
 224     //splitting along only dimension
 225     _local_array_indices.resize(1);
 226     _nb_procs_per_dim.resize(1);
 227
 228     _local_array_indices[0].resize(_proc_group->size()+1);
 229     _local_array_indices[0][0]=0;
 230     _nb_procs_per_dim[0]=_proc_group->size();
 231
 232     for (int i=1; i<=_proc_group->size(); i++)
 233       {
 234         _local_array_indices[0][i]=_local_array_indices[0][i-1]+
 235           nbelems_per_proc[i-1];
 236         _nb_elems+=nbelems_per_proc[i-1];
 237       }
 238     _cycle_type.resize(1);
 239     _cycle_type[0]=ParaMEDMEM::Block;
 240     delete[] nbelems_per_proc;
 241   }
 242
 243   BlockTopology::~BlockTopology()
 244   {
 245     if (_owns_processor_group)
 246       delete _proc_group;
 247   }
 248
 249   /*! Retrieves the min and max indices of the domain stored locally
 250    * for each dimension. The output vector has the topology dimension
 251    * as a size and each pair <int,int> contains min and max. Indices
 252    * range from min to max-1.
 253    */
 254   std::vector<std::pair<int,int> > BlockTopology::getLocalArrayMinMax() const
 255   {
 256     vector<pair<int,int> > local_indices (_dimension);
 257     int myrank=_proc_group->myRank();
 258     int increment=1;
 259     for (int i=_dimension-1; i>=0; i--)
 260       {
 261         increment *=_nb_procs_per_dim[i];
 262         int idim=myrank%increment;
 263         local_indices[i].first=_local_array_indices[i][idim];
 264         local_indices[i].second=_local_array_indices[i][idim+1];
 265         cout << local_indices[i].first << " "<< local_indices[i].second<<endl;
 266       }
 267     return local_indices;
 268   }
 269
 270   /*! Serializes the data contained in the Block Topology
 271    * for communication purposes*/
 272   void BlockTopology::serialize(int* & serializer, int& size) const
 273   {
 274     vector<int> buffer;
 275
 276     buffer.push_back(_dimension);
 277     buffer.push_back(_nb_elems);
 278     for (int i=0; i<_dimension; i++)
 279       {
 280         buffer.push_back(_nb_procs_per_dim[i]);
 281         buffer.push_back(_cycle_type[i]);
 282         buffer.push_back(_local_array_indices[i].size());
 283         for (int j=0; j<(int)_local_array_indices[i].size(); j++)
 284           buffer.push_back(_local_array_indices[i][j]);
 285       }
 286
 287     //serializing the comm group
 288     int size_comm=_proc_group->size();
 289     buffer.push_back(size_comm);
 290     MPIProcessorGroup world_group(_proc_group->getCommInterface());
 291     for (int i=0; i<size_comm;i++)
 292       {
 293         int world_rank=world_group.translateRank(_proc_group, i);
 294         buffer.push_back(world_rank);
 295       }
 296
 297     serializer=new int[buffer.size()];
 298     size=buffer.size();
 299     copy(buffer.begin(), buffer.end(), serializer);
 300   }
 301
 302   /*!
 303    *
 304    * Unserializes the data contained in the Block Topology
 305    * after communication. Uses the same structure as the one used for serialize()
 306    *
 307    */
 308   void BlockTopology::unserialize(const int* serializer,const CommInterface& comm_interface)
 309   {
 310     const int* ptr_serializer=serializer;
 311     cout << "unserialize..."<<endl;
 312     _dimension=*(ptr_serializer++);
 313     cout << "dimension "<<_dimension<<endl;
 314     _nb_elems=*(ptr_serializer++);
 315     cout << "nbelems "<<_nb_elems<<endl;
 316     _nb_procs_per_dim.resize(_dimension);
 317     _cycle_type.resize(_dimension);
 318     _local_array_indices.resize(_dimension);
 319     for (int i=0; i<_dimension; i++)
 320       {
 321         _nb_procs_per_dim[i]=*(ptr_serializer++);
 322         _cycle_type[i]=(CYCLE_TYPE)*(ptr_serializer++);
 323         _local_array_indices[i].resize(*(ptr_serializer++));
 324         for (int j=0; j<(int)_local_array_indices[i].size(); j++)
 325           _local_array_indices[i][j]=*(ptr_serializer++);
 326       }
 327     set<int> procs;
 328     int size_comm=*(ptr_serializer++);
 329     for (int i=0; i<size_comm; i++)
 330       procs.insert(*(ptr_serializer++));
 331     cout << "unserialize..."<<procs.size()<<endl;
 332     _proc_group=new MPIProcessorGroup(comm_interface,procs);
 333     _owns_processor_group=true;
 334     //TODO manage memory ownership of _proc_group
 335   }
 336 }