src/ParaMEDMEM/DisjointDEC.cxx

   1 // Copyright (C) 2007-2023  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 "DisjointDEC.hxx"
  21 #include "CommInterface.hxx"
  22 #include "Topology.hxx"
  23 #include "BlockTopology.hxx"
  24 #include "ComponentTopology.hxx"
  25 #include "ParaFIELD.hxx"
  26 #include "ParaMESH.hxx"
  27 #include "ICoCoField.hxx"
  28 #include "ICoCoMEDDoubleField.hxx"
  29 #include "MPIProcessorGroup.hxx"
  30
  31 #include <cmath>
  32 #include <iostream>
  33
  34
  35 namespace MEDCoupling
  36 {
  37   DisjointDEC::DisjointDEC(ProcessorGroup& source_group, ProcessorGroup& target_group):
  38       _local_field(0),
  39       _source_group(&source_group),
  40       _target_group(&target_group),
  41       _comm_interface(0),
  42       _owns_field(false),
  43       _owns_groups(false),
  44       _union_comm(MPI_COMM_NULL)
  45   {
  46     checkPartitionGroup();
  47     _union_group = source_group.fuse(target_group);
  48   }
  49
  50   DisjointDEC::DisjointDEC(const DisjointDEC& s):
  51       DEC(s),
  52       _local_field(0),
  53       _union_group(0),
  54       _source_group(0),
  55       _target_group(0),
  56       _comm_interface(0),
  57       _owns_field(false),
  58       _owns_groups(false),
  59       _union_comm(MPI_COMM_NULL)
  60   {
  61     copyInstance(s);
  62   }
  63
  64   DisjointDEC & DisjointDEC::operator=(const DisjointDEC& s)
  65   {
  66     cleanInstance();
  67     copyInstance(s);
  68     return *this;
  69
  70   }
  71
  72   void DisjointDEC::copyInstance(const DisjointDEC& other)
  73   {
  74     DEC::copyFrom(other);
  75     if (other._union_comm != MPI_COMM_NULL)
  76       {
  77         // Tricky: the DEC is responsible for the management of _union_comm. And this comm is referenced by
  78         // the MPIProcGroups (source/targets). In the case where _union_comm is not NULL we must take care of rebuilding the
  79         // MPIProcGroups with a communicator that will survive the destruction of 'other'.
  80         _owns_groups = true;
  81         MPI_Comm_dup(other._union_comm, &_union_comm);
  82 //        std::cout << "DUP union comm - new is "<< _union_comm << "\n";
  83         _target_group = new MPIProcessorGroup(*_comm_interface, other._target_group->getProcIDs(), _union_comm);
  84         _source_group = new MPIProcessorGroup(*_comm_interface, other._source_group->getProcIDs(), _union_comm);
  85       }
  86     else{
  87       if(other._target_group)
  88         {
  89           _target_group=other._target_group->deepCopy();
  90           _owns_groups=true;
  91         }
  92       if(other._source_group)
  93         {
  94           _source_group=other._source_group->deepCopy();
  95           _owns_groups=true;
  96         }
  97     }
  98     if (_source_group && _target_group)
  99       _union_group = _source_group->fuse(*_target_group);
 100   }
 101
 102   DisjointDEC::DisjointDEC(const std::set<int>& source_ids,
 103                            const std::set<int>& target_ids,
 104                            const MPI_Comm& world_comm):
 105      _local_field(0),
 106      _comm_interface(0),
 107      _owns_field(false),
 108      _owns_groups(true),
 109      _union_comm(MPI_COMM_NULL)
 110   {
 111     MEDCoupling::CommInterface comm;
 112     // Create the list of procs including source and target
 113     std::set<int> union_ids; // source and target ids in world_comm
 114     union_ids.insert(source_ids.begin(),source_ids.end());
 115     union_ids.insert(target_ids.begin(),target_ids.end());
 116     if(union_ids.size()!=(source_ids.size()+target_ids.size()))
 117       throw INTERP_KERNEL::Exception("DisjointDEC constructor : source_ids and target_ids overlap partially or fully. This type of DEC does not support it! OverlapDEC class could be the solution!");
 118     int* union_ranks_world=new int[union_ids.size()]; // ranks of sources and targets in world_comm
 119     std::copy(union_ids.begin(), union_ids.end(), union_ranks_world);
 120
 121     // Create a communicator on these procs
 122     MPI_Group union_group,world_group;
 123     comm.commGroup(world_comm,&world_group);
 124     comm.groupIncl(world_group,(int)union_ids.size(),union_ranks_world,&union_group);
 125     comm.commCreate(world_comm,union_group,&_union_comm);
 126     delete[] union_ranks_world;
 127     if (_union_comm==MPI_COMM_NULL)
 128       { // This process is not in union
 129         _source_group=0;
 130         _target_group=0;
 131         _union_group=0;
 132         comm.groupFree(&union_group);
 133         comm.groupFree(&world_group);
 134         return;
 135       }
 136
 137     // Translate source_ids and target_ids from world_comm to union_comm
 138     int* source_ranks_world=new int[source_ids.size()]; // ranks of sources in world_comm
 139     std::copy(source_ids.begin(), source_ids.end(),source_ranks_world);
 140     int* source_ranks_union=new int[source_ids.size()]; // ranks of sources in union_comm
 141     int* target_ranks_world=new int[target_ids.size()]; // ranks of targets in world_comm
 142     std::copy(target_ids.begin(), target_ids.end(),target_ranks_world);
 143     int* target_ranks_union=new int[target_ids.size()]; // ranks of targets in union_comm
 144     MPI_Group_translate_ranks(world_group,(int)source_ids.size(),source_ranks_world,union_group,source_ranks_union);
 145     MPI_Group_translate_ranks(world_group,(int)target_ids.size(),target_ranks_world,union_group,target_ranks_union);
 146     std::set<int> source_ids_union;
 147     for (int i=0;i<(int)source_ids.size();i++)
 148       source_ids_union.insert(source_ranks_union[i]);
 149     std::set<int> target_ids_union;
 150     for (int i=0;i<(int)target_ids.size();i++)
 151       target_ids_union.insert(target_ranks_union[i]);
 152     delete [] source_ranks_world;
 153     delete [] source_ranks_union;
 154     delete [] target_ranks_world;
 155     delete [] target_ranks_union;
 156
 157     // Create the MPIProcessorGroups
 158     _source_group = new MPIProcessorGroup(comm,source_ids_union,_union_comm);
 159     _target_group = new MPIProcessorGroup(comm,target_ids_union,_union_comm);
 160     _union_group = _source_group->fuse(*_target_group);
 161     comm.groupFree(&union_group);
 162     comm.groupFree(&world_group);
 163   }
 164
 165   DisjointDEC::~DisjointDEC()
 166   {
 167     cleanInstance();
 168   }
 169
 170   void DisjointDEC::cleanInstance()
 171   {
 172     if(_owns_field)
 173       {
 174         delete _local_field;
 175       }
 176     _local_field=nullptr;
 177     _owns_field=false;
 178     if(_owns_groups)
 179       {
 180         delete _source_group;
 181         delete _target_group;
 182       }
 183     _owns_groups=false;
 184     _source_group=nullptr;
 185     _target_group=nullptr;
 186     delete _union_group;
 187     _union_group=nullptr;
 188     if (_union_comm != MPI_COMM_NULL)
 189       _comm_interface->commFree(&_union_comm);
 190     _union_comm = MPI_COMM_NULL;
 191   }
 192
 193   /**
 194    * Check that the sources and targets procs form a partition of the world communicator referenced in the groups.
 195    * This world communicator is not necessarily MPI_WORLD_COMM, but it has to be covered completely for the DECs to work.
 196    */
 197   void DisjointDEC::checkPartitionGroup() const
 198   {
 199     int size = -1;
 200     MPIProcessorGroup * tgt = static_cast<MPIProcessorGroup *>(_target_group);
 201     MPIProcessorGroup * src = static_cast<MPIProcessorGroup *>(_source_group);
 202     MPI_Comm comm_t = tgt->getWorldComm();
 203     MPI_Comm comm_s = src->getWorldComm();
 204     if (comm_t != comm_s)
 205       throw INTERP_KERNEL::Exception("DisjointDEC constructor: Inconsistent world communicator when building DisjointDEC");
 206     MPI_Comm_size(comm_t, &size);
 207
 208     std::set<int> union_ids; // source and target ids in world_comm
 209     union_ids.insert(src->getProcIDs().begin(),src->getProcIDs().end());
 210     union_ids.insert(tgt->getProcIDs().begin(),tgt->getProcIDs().end());
 211     if((int)union_ids.size()!=size)
 212       throw INTERP_KERNEL::Exception("DisjointDEC constructor: source_ids and target_ids do not form a partition of the communicator! Restrain the world communicator passed to MPIProcessorGroup ctor.");
 213   }
 214
 215   void DisjointDEC::setNature(NatureOfField nature)
 216   {
 217     if(_local_field)
 218       _local_field->getField()->setNature(nature);
 219   }
 220
 221   /*! Attaches a local field to a DEC.
 222     If the processor is on the receiving end of the DEC, the field
 223     will be updated by a recvData() call.
 224     Reversely, if the processor is on the sending end, the field will be read, possibly transformed, and sent appropriately to the other side.
 225   */
 226   void DisjointDEC::attachLocalField(const ParaFIELD *field, bool ownPt)
 227   {
 228     if(!isInUnion())
 229       return ;
 230     if(_owns_field)
 231       delete _local_field;
 232     _local_field=field;
 233     _owns_field=ownPt;
 234     _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface());
 235     compareFieldAndMethod();
 236   }
 237
 238   /*! Attaches a local field to a DEC. The method will test whether the processor
 239     is on the source or the target side and will associate the mesh underlying the
 240     field to the local side.
 241
 242     If the processor is on the receiving end of the DEC, the field
 243     will be updated by a recvData() call.
 244     Reversely, if the processor is on the sending end, the field will be read, possibly transformed,
 245     and sent appropriately to the other side.
 246   */
 247
 248   void DisjointDEC::attachLocalField(MEDCouplingFieldDouble *field)
 249   {
 250     if(!isInUnion())
 251       return ;
 252     ProcessorGroup* local_group;
 253     if (_source_group->containsMyRank())
 254       local_group=_source_group;
 255     else if (_target_group->containsMyRank())
 256       local_group=_target_group;
 257     else
 258       throw INTERP_KERNEL::Exception("Invalid procgroup for field attachment to DEC");
 259     ParaMESH *paramesh=new ParaMESH(static_cast<MEDCouplingPointSet *>(const_cast<MEDCouplingMesh *>(field->getMesh())),*local_group,field->getMesh()->getName());
 260     ParaFIELD *tmp=new ParaFIELD(field, paramesh, *local_group);
 261     tmp->setOwnSupport(true);
 262     attachLocalField(tmp,true);
 263     //_comm_interface=&(local_group->getCommInterface());
 264   }
 265
 266   /*!
 267     Attaches a local field to a DEC.
 268     If the processor is on the receiving end of the DEC, the field
 269     will be updated by a recvData() call.
 270     Reversely, if the processor is on the sending end, the field will be read, possibly transformed, and sent appropriately to the other side.
 271     The field type is a generic ICoCo Field, so that the DEC can couple a number of different fields :
 272     - a ICoCo::MEDDoubleField, that is created from a MEDCoupling structure
 273
 274   */
 275   void DisjointDEC::attachLocalField(const ICoCo::MEDDoubleField *field)
 276   {
 277     if(!isInUnion())
 278       return ;
 279     if(!field)
 280       throw INTERP_KERNEL::Exception("DisjointDEC::attachLocalField : ICoCo::MEDDoubleField pointer is NULL !");
 281     attachLocalField(field->getMCField());
 282   }
 283
 284   /*!
 285     Computes the field norm over its support
 286     on the source side and renormalizes the field on the target side
 287     so that the norms match.
 288
 289     \f[
 290     I_{source}=\sum_{i=1}^{n_{source}}V_{i}.|\Phi^{source}_{i}|^2,
 291     \f]
 292
 293     \f[
 294     I_{target}=\sum_{i=1}^{n_{target}}V_{i}.|\Phi^{target}_{i}|^2,
 295     \f]
 296
 297     \f[
 298     \Phi^{target}:=\Phi^{target}.\sqrt{I_{source}/I_{target}}.
 299     \f]
 300
 301   */
 302   void DisjointDEC::renormalizeTargetField(bool isWAbs)
 303   {
 304     if (_source_group->containsMyRank())
 305       for (int icomp=0; icomp<(int)_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
 306         {
 307           double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
 308           double source_norm = total_norm;
 309           _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
 310
 311         }
 312     if (_target_group->containsMyRank())
 313       {
 314         for (int icomp=0; icomp<(int)_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
 315           {
 316             double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
 317             double source_norm=total_norm;
 318             _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
 319
 320             if (fabs(total_norm)>1e-100)
 321               _local_field->getField()->applyLin(source_norm/total_norm,0.0,icomp+1);
 322           }
 323       }
 324   }
 325
 326   bool DisjointDEC::isInSourceSide() const
 327   {
 328     if(!_source_group)
 329       return false;
 330     return _source_group->containsMyRank();
 331   }
 332
 333   bool DisjointDEC::isInTargetSide() const
 334   {
 335     if(!_target_group)
 336       return false;
 337     return _target_group->containsMyRank();
 338   }
 339
 340   bool DisjointDEC::isInUnion() const
 341   {
 342     if(!_union_group)
 343       return false;
 344     return _union_group->containsMyRank();
 345   }
 346
 347   void DisjointDEC::compareFieldAndMethod() const
 348   {
 349     if (_local_field)
 350       {
 351         TypeOfField entity = _local_field->getField()->getTypeOfField();
 352         if ( getMethod() == "P0" )
 353           {
 354             if ( entity != ON_CELLS )
 355               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
 356                                              " For P0 interpolation, field must be on MED_CELL's");
 357           }
 358         else if ( getMethod() == "P1" )
 359           {
 360             if ( entity != ON_NODES )
 361               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
 362                                              " For P1 interpolation, field must be on MED_NODE's");
 363           }
 364         else if ( getMethod() == "P1d" )
 365           {
 366             if ( entity != ON_CELLS )
 367               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
 368                                              " For P1d interpolation, field must be on MED_CELL's");
 369             if ( _target_group->containsMyRank() )
 370               throw INTERP_KERNEL::Exception("Projection to P1d field not supported");
 371           }
 372         else
 373           {
 374             throw INTERP_KERNEL::Exception("Unknown interpolation method. Possible methods: P0, P1, P1d");
 375           }
 376       }
 377   }
 378
 379   /*!
 380     If way==true, source procs call sendData() and target procs call recvData().
 381     if way==false, it's the other way round.
 382   */
 383   void DisjointDEC::sendRecvData(bool way)
 384   {
 385     if(!isInUnion())
 386       return;
 387     if(isInSourceSide())
 388       {
 389         if(way)
 390           sendData();
 391         else
 392           recvData();
 393       }
 394     else if(isInTargetSide())
 395       {
 396         if(way)
 397           recvData();
 398         else
 399           sendData();
 400       }
 401   }
 402 }