src/ParaMEDMEM/BlockTopology.cxx

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