src/ParaMEDMEM/DisjointDEC.cxx

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