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