1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #include "MEDPARTITIONER_ParaDomainSelector.hxx"
21 #include "MEDPARTITIONER_UserGraph.hxx"
22 #include "MEDPARTITIONER_Utils.hxx"
24 #include "MEDCouplingUMesh.hxx"
25 #include "MEDCouplingSkyLineArray.hxx"
35 * \brief Constructor. Find out my rank and world size
37 MEDPARTITIONER::ParaDomainSelector::ParaDomainSelector(bool mesure_memory)
38 :_rank(0),_world_size(1), _nb_result_domains(-1), _init_time(0.0),
39 _mesure_memory(mesure_memory), _init_memory(0), _max_memory(0)
42 if (MyGlobals::_Rank==-1)
44 MPI_Init(0,0); //do once only
45 MPI_Comm_size(MPI_COMM_WORLD,&_world_size) ;
46 MPI_Comm_rank(MPI_COMM_WORLD,&_rank) ;
50 _world_size=MyGlobals::_World_Size;
51 _rank=MyGlobals::_Rank;
53 _init_time = MPI_Wtime();
58 if (MyGlobals::_Verbose>10)
59 std::cout << "WARNING : ParaDomainSelector contructor without parallel_mode World_Size=1 by default" << std::endl;
61 MyGlobals::_World_Size=_world_size;
62 MyGlobals::_Rank=_rank;
64 if (MyGlobals::_Verbose>200) std::cout << "proc " << MyGlobals::_Rank << " of " << MyGlobals::_World_Size << std::endl;
68 MEDPARTITIONER::ParaDomainSelector::~ParaDomainSelector()
73 * \brief Return true if is running on different hosts
75 bool MEDPARTITIONER::ParaDomainSelector::isOnDifferentHosts() const
78 if ( _world_size < 2 ) return false;
81 char name_here[ MPI_MAX_PROCESSOR_NAME+1 ], name_there[ MPI_MAX_PROCESSOR_NAME+1 ];
83 MPI_Get_processor_name( name_here, &size);
85 int next_proc = (rank() + 1) % nbProcs();
86 int prev_proc = (rank() - 1 + nbProcs()) % nbProcs();
90 MPI_Sendrecv((void*)&name_here[0], MPI_MAX_PROCESSOR_NAME, MPI_CHAR, next_proc, tag,
91 (void*)&name_there[0], MPI_MAX_PROCESSOR_NAME, MPI_CHAR, prev_proc, tag,
92 MPI_COMM_WORLD, &status);
94 //bug: (isOnDifferentHosts here and there) is not (isOnDifferentHosts somewhere)
95 //return string(name_here) != string(name_there);
99 if (std::string(name_here) != std::string(name_there))
101 MPI_Allreduce( &same, &sum_same, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );
102 return (sum_same != nbProcs());
108 * \brief Return true if the domain with domainIndex is to be loaded on this proc
109 * \param domainIndex - index of mesh domain
110 * \retval bool - to load or not
112 bool MEDPARTITIONER::ParaDomainSelector::isMyDomain(int domainIndex) const
115 return (_rank == getProcessorID( domainIndex ));
119 * \brief Return processor id where the domain with domainIndex resides
120 * \param domainIndex - index of mesh domain
121 * \retval int - processor id
123 int MEDPARTITIONER::ParaDomainSelector::getProcessorID(int domainIndex) const
126 return ( domainIndex % _world_size );
130 * \brief Gather info on nb of cell entities on each processor and return total nb.
133 * 1) for MED_CELL to know global id shift for domains at graph construction;
134 * 2) for sub-entity to know total nb of sub-entities before creating those of joints
136 void MEDPARTITIONER::ParaDomainSelector::gatherNbOf(const std::vector<MEDCoupling::MEDCouplingUMesh*>& domain_meshes)
139 // get nb of elems of each domain mesh
140 int nb_domains=domain_meshes.size();
141 std::vector<int> nb_elems(nb_domains*2, 0); //NumberOfCells & NumberOfNodes
142 for (int i=0; i<nb_domains; ++i)
143 if ( domain_meshes[i] )
145 nb_elems[i*2] = domain_meshes[i]->getNumberOfCells();
146 nb_elems[i*2+1] = domain_meshes[i]->getNumberOfNodes();
148 // receive nb of elems from other procs
149 std::vector<int> all_nb_elems;
150 if (MyGlobals::_World_Size==1)
152 all_nb_elems=nb_elems;
157 all_nb_elems.resize( nb_domains*2 );
158 MPI_Allreduce((void*)&nb_elems[0], (void*)&all_nb_elems[0], nb_domains*2, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
160 throw INTERP_KERNEL::Exception("not(HAVE_MPI) incompatible with MPI_World_Size>1");
163 int total_nb_cells=0, total_nb_nodes=0;
164 for (int i=0; i<nb_domains; ++i)
166 total_nb_cells+=all_nb_elems[i*2];
167 total_nb_nodes+=all_nb_elems[i*2+1];
170 if (MyGlobals::_Is0verbose>10)
171 std::cout << "totalNbCells " << total_nb_cells << " totalNbNodes " << total_nb_nodes << std::endl;
173 std::vector<int>& cell_shift_by_domain=_cell_shift_by_domain;
174 std::vector<int>& node_shift_by_domain=_node_shift_by_domain;
175 std::vector<int>& face_shift_by_domain=_face_shift_by_domain;
177 std::vector< int > ordered_nbs_cell, ordered_nbs_node, domain_order(nb_domains);
178 ordered_nbs_cell.push_back(0);
179 ordered_nbs_node.push_back(0);
180 for (int iproc=0; iproc<nbProcs(); ++iproc)
181 for (int idomain=0; idomain<nb_domains; ++idomain)
182 if (getProcessorID( idomain )==iproc)
184 domain_order[idomain] = ordered_nbs_cell.size() - 1;
185 ordered_nbs_cell.push_back( ordered_nbs_cell.back() + all_nb_elems[idomain*2] );
186 ordered_nbs_node.push_back( ordered_nbs_node.back() + all_nb_elems[idomain*2+1] );
188 cell_shift_by_domain.resize( nb_domains+1, 0 );
189 node_shift_by_domain.resize( nb_domains+1, 0 );
190 face_shift_by_domain.resize( nb_domains+1, 0 );
191 for (int idomain=0; idomain<nb_domains; ++idomain)
193 cell_shift_by_domain[ idomain ] = ordered_nbs_cell[ domain_order[ idomain ]];
194 node_shift_by_domain[ idomain ] = ordered_nbs_node[ domain_order[ idomain ]];
196 cell_shift_by_domain.back() = ordered_nbs_cell.back(); // to know total nb of elements
197 node_shift_by_domain.back() = ordered_nbs_node.back(); // to know total nb of elements
199 if (MyGlobals::_Is0verbose>300)
201 std::cout << "proc " << MyGlobals::_Rank << " : cellShiftByDomain ";
202 for (int i=0; i<=nb_domains; ++i)
203 std::cout << cell_shift_by_domain[i] << "|";
204 std::cout << std::endl;
205 std::cout << "proc " << MyGlobals::_Rank << " : nodeShiftBy_domain ";
206 for (int i=0; i<=nb_domains; ++i)
207 std::cout << node_shift_by_domain[i] << "|";
208 std::cout << std::endl;
210 // fill _nb_vert_of_procs (is Vtxdist)
211 _nb_vert_of_procs.resize(_world_size+1);
212 _nb_vert_of_procs[0] = 0; // base = 0
213 for (int i=0; i<nb_domains; ++i)
215 int rankk = getProcessorID(i);
216 _nb_vert_of_procs[rankk+1] += all_nb_elems[i*2];
218 for (std::size_t i=1; i<_nb_vert_of_procs.size(); ++i)
219 _nb_vert_of_procs[i] += _nb_vert_of_procs[i-1]; // to CSR format : cumulated
221 if (MyGlobals::_Is0verbose>200)
223 std::cout << "proc " << MyGlobals::_Rank << " : gatherNbOf : vtxdist is ";
224 for (int i = 0; i <= _world_size; ++i)
225 std::cout << _nb_vert_of_procs[i] << " ";
226 std::cout << std::endl;
234 * \brief Return distribution of the graph vertices among the processors
235 * \retval int* - array containing nb of vertices (=cells) on all processors
237 * gatherNbOf() must be called before.
238 * The result array is to be used as the first arg of ParMETIS_V3_PartKway() and
239 * is freed by ParaDomainSelector.
241 int *MEDPARTITIONER::ParaDomainSelector::getProcVtxdist() const
244 if (_nb_vert_of_procs.empty())
245 throw INTERP_KERNEL::Exception("_nb_vert_of_procs not set");
246 return const_cast<int*>(& _nb_vert_of_procs[0]);
250 * \brief Return nb of cells in domains with lower index.
252 * gatherNbOf() must be called before.
253 * Result added to local id on given domain gives id in the whole distributed mesh
255 int MEDPARTITIONER::ParaDomainSelector::getDomainCellShift(int domainIndex) const
258 if (_cell_shift_by_domain.empty())
259 throw INTERP_KERNEL::Exception("_cell_shift_by_domain not set");
260 return _cell_shift_by_domain[domainIndex];
263 int MEDPARTITIONER::ParaDomainSelector::getDomainNodeShift(int domainIndex) const
266 if (_node_shift_by_domain.empty())
267 throw INTERP_KERNEL::Exception("_node_shift_by_domain not set");
268 return _node_shift_by_domain[domainIndex];
272 * \brief Return nb of nodes on processors with lower rank.
274 * gatherNbOf() must be called before.
275 * Result added to global id on this processor gives id in the whole distributed mesh
277 int MEDPARTITIONER::ParaDomainSelector::getProcNodeShift() const
280 if (_nb_vert_of_procs.empty())
281 throw INTERP_KERNEL::Exception("_nb_vert_of_procs not set");
282 return _nb_vert_of_procs[_rank];
286 * \brief Gather graphs from all processors into one
288 std::auto_ptr<MEDPARTITIONER::Graph> MEDPARTITIONER::ParaDomainSelector::gatherGraph(const Graph* graph) const
290 Graph* glob_graph = 0;
299 std::vector<int> index_size_of_proc( nbProcs() ); // index sizes - 1
300 for ( std::size_t i = 1; i < _nb_vert_of_procs.size(); ++i )
301 index_size_of_proc[i-1] = _nb_vert_of_procs[ i ] - _nb_vert_of_procs[ i-1 ];
303 int index_size = 1 + _cell_shift_by_domain.back();
304 int *graph_index = new int[ index_size ];
305 const int *index = graph->getGraph()->getIndex();
306 int *proc_index_displacement = const_cast<int*>( & _nb_vert_of_procs[0] );
308 MPI_Allgatherv((void*) (index+1), // send local index except first 0 (or 1)
309 index_size_of_proc[_rank], // index size on this proc
311 (void*) graph_index, // receive indices
312 & index_size_of_proc[0], // index size on each proc
313 proc_index_displacement, // displacement of each proc index
316 graph_index[0] = index[0]; // it is not overwritten thanks to proc_index_displacement[0]==1
318 // get sizes of graph values on each proc by the got indices of graphs
319 std::vector< int > value_size_of_proc( nbProcs() ), proc_value_displacement(1,0);
320 for ( int i = 0; i < nbProcs(); ++i )
322 if ( index_size_of_proc[i] > 0 )
323 value_size_of_proc[i] = graph_index[ proc_index_displacement[ i+1 ]-1 ] - graph_index[0];
325 value_size_of_proc[i] = 0;
326 proc_value_displacement.push_back( proc_value_displacement.back() + value_size_of_proc[i] );
329 // update graph_index
330 for ( int i = 1; i < nbProcs(); ++i )
332 int shift = graph_index[ proc_index_displacement[i]-1 ]-graph_index[0];
333 for ( int j = proc_index_displacement[i]; j < proc_index_displacement[i+1]; ++j )
334 graph_index[ j ] += shift;
341 int value_size = graph_index[ index_size-1 ] - graph_index[ 0 ];
342 int *graph_value = new int[ value_size ];
343 const int *value = graph->getGraph()->getValues();
345 MPI_Allgatherv((void*) value, // send local value
346 value_size_of_proc[_rank], // value size on this proc
348 (void*) graph_value, // receive values
349 & value_size_of_proc[0], // value size on each proc
350 & proc_value_displacement[0], // displacement of each proc value
358 int * partition = new int[ _cell_shift_by_domain.back() ];
359 const int* part = graph->getPart();
361 MPI_Allgatherv((void*) part, // send local partition
362 index_size_of_proc[_rank], // index size on this proc
364 (void*)(partition-1), // -1 compensates proc_index_displacement[0]==1
365 & index_size_of_proc[0], // index size on each proc
366 proc_index_displacement, // displacement of each proc index
374 // MEDCouplingSkyLineArray* array =
375 // new MEDCouplingSkyLineArray( index_size-1, value_size, graph_index, graph_value, true );
377 // glob_graph = new UserGraph( array, partition, index_size-1 );
385 return std::auto_ptr<Graph>( glob_graph );
390 * \brief Set nb of cell/cell pairs in a joint between domains
392 void MEDPARTITIONER::ParaDomainSelector::setNbCellPairs( int nb_cell_pairs, int dist_domain, int loc_domain )
394 // This method is needed for further computing global numbers of faces in joint.
395 // Store if both domains are on this proc else on one of procs only
396 if ( isMyDomain( dist_domain ) || dist_domain < loc_domain )
398 if ( _nb_cell_pairs_by_joint.empty() )
399 _nb_cell_pairs_by_joint.resize( _nb_result_domains*(_nb_result_domains+1), 0);
401 int joint_id = jointId( loc_domain, dist_domain );
402 _nb_cell_pairs_by_joint[ joint_id ] = nb_cell_pairs;
407 //================================================================================
409 * \brief Return nb of cell/cell pairs in a joint between domains on different procs
411 //================================================================================
413 int MEDPARTITIONER::ParaDomainSelector::getNbCellPairs( int dist_domain, int loc_domain ) const
417 int joint_id = jointId( loc_domain, dist_domain );
418 return _nb_cell_pairs_by_joint[ joint_id ];
421 //================================================================================
423 * \brief Gather size of each joint
425 //================================================================================
427 void MEDPARTITIONER::ParaDomainSelector::gatherNbCellPairs()
429 if ( _nb_cell_pairs_by_joint.empty() )
430 _nb_cell_pairs_by_joint.resize( _nb_result_domains*(_nb_result_domains+1), 0);
433 std::vector< int > send_buf = _nb_cell_pairs_by_joint;
435 MPI_Allreduce((void*)&send_buf[0],
436 (void*)&_nb_cell_pairs_by_joint[0],
437 _nb_cell_pairs_by_joint.size(),
438 MPI_INT, MPI_SUM, MPI_COMM_WORLD);
440 // check that the set nbs of cell pairs are correct,
441 // namely that each joint is treated on one proc only
442 for ( std::size_t j = 0; j < _nb_cell_pairs_by_joint.size(); ++j )
443 if ( _nb_cell_pairs_by_joint[j] != send_buf[j] && send_buf[j]>0 )
444 throw INTERP_KERNEL::Exception("invalid nb of cell pairs");
447 //================================================================================
449 * \brief Return the first global id of sub-entity for the joint
451 //================================================================================
453 int MEDPARTITIONER::ParaDomainSelector::getFisrtGlobalIdOfSubentity( int loc_domain, int dist_domain ) const
455 // total_nb_faces includes faces existing before creation of joint faces
456 // (got in gatherNbOf( MED_FACE )).
459 int total_nb_faces = _face_shift_by_domain.empty() ? 0 : _face_shift_by_domain.back();
460 int id = total_nb_faces + 1;
462 if ( _nb_cell_pairs_by_joint.empty() )
463 throw INTERP_KERNEL::Exception("gatherNbCellPairs() must be called before");
464 int joint_id = jointId( loc_domain, dist_domain );
465 for ( int j = 0; j < joint_id; ++j )
466 id += _nb_cell_pairs_by_joint[ j ];
471 //================================================================================
473 * \brief Send-receive local ids of joint faces
475 //================================================================================
477 int *MEDPARTITIONER::ParaDomainSelector::exchangeSubentityIds( int loc_domain, int dist_domain,
478 const std::vector<int>& loc_ids_here ) const
480 int* loc_ids_dist = new int[ loc_ids_here.size()];
482 int dest = getProcessorID( dist_domain );
483 int tag = 2002 + jointId( loc_domain, dist_domain );
485 MPI_Sendrecv((void*)&loc_ids_here[0], loc_ids_here.size(), MPI_INT, dest, tag,
486 (void*) loc_ids_dist, loc_ids_here.size(), MPI_INT, dest, tag,
487 MPI_COMM_WORLD, &status);
494 //================================================================================
496 * \brief Return identifier for a joint
498 //================================================================================
500 int MEDPARTITIONER::ParaDomainSelector::jointId( int local_domain, int distant_domain ) const
503 if (_nb_result_domains < 0)
504 throw INTERP_KERNEL::Exception("setNbDomains() must be called before");
506 if ( local_domain < distant_domain )
507 std::swap( local_domain, distant_domain );
508 return local_domain * _nb_result_domains + distant_domain;
512 //================================================================================
514 * \brief Return time passed from construction in seconds
516 //================================================================================
518 double MEDPARTITIONER::ParaDomainSelector::getPassedTime() const
521 return MPI_Wtime() - _init_time;
528 Sends content of \a mesh to processor \a target. To be used with \a recvMesh method.
529 \param mesh mesh to be sent
530 \param target processor id of the target
533 void MEDPARTITIONER::ParaDomainSelector::sendMesh(const MEDCoupling::MEDCouplingUMesh& mesh, int target) const
536 throw INTERP_KERNEL::Exception("ParaDomainSelector::sendMesh : incoherent call in non_MPI mode");
538 if (MyGlobals::_Verbose>600)
539 std::cout << "proc " << _rank << " : sendMesh '" << mesh.getName() << "' size " << mesh.getNumberOfCells() << " to " << target << std::endl;
540 // First stage : sending sizes
541 // ------------------------------
542 std::vector<int> tinyInfoLocal;
543 std::vector<std::string> tinyInfoLocalS;
544 std::vector<double> tinyInfoLocalD;
545 //Getting tiny info of local mesh to allow the distant proc to initialize and allocate
546 //the transmitted mesh.
547 mesh.getTinySerializationInformation(tinyInfoLocalD,tinyInfoLocal,tinyInfoLocalS);
548 tinyInfoLocal.push_back(mesh.getNumberOfCells());
549 int tinySize=tinyInfoLocal.size();
550 MPI_Send(&tinySize, 1, MPI_INT, target, 1113, MPI_COMM_WORLD);
551 MPI_Send(&tinyInfoLocal[0], tinyInfoLocal.size(), MPI_INT, target, 1112, MPI_COMM_WORLD);
553 if (mesh.getNumberOfCells()>0) //no sends if empty
555 MEDCoupling::DataArrayInt *v1Local=0;
556 MEDCoupling::DataArrayDouble *v2Local=0;
557 //serialization of local mesh to send data to distant proc.
558 mesh.serialize(v1Local,v2Local);
561 if(v1Local) //if empty getNbOfElems() is 1!
563 nbLocalElems=v1Local->getNbOfElems(); // if empty be 1!
564 ptLocal=v1Local->getPointer();
566 MPI_Send(ptLocal, nbLocalElems, MPI_INT, target, 1111, MPI_COMM_WORLD);
569 if(v2Local) //if empty be 0!
571 nbLocalElems2=v2Local->getNbOfElems();
572 ptLocal2=v2Local->getPointer();
574 MPI_Send(ptLocal2, nbLocalElems2, MPI_DOUBLE, target, 1110, MPI_COMM_WORLD);
575 if(v1Local) v1Local->decrRef();
576 if(v2Local) v2Local->decrRef();
581 /*! Receives messages from proc \a source to fill mesh \a mesh.
582 To be used with \a sendMesh method.
583 \param mesh pointer to mesh that is filled
584 \param source processor id of the incoming messages
586 void MEDPARTITIONER::ParaDomainSelector::recvMesh(MEDCoupling::MEDCouplingUMesh*& mesh, int source)const
589 throw INTERP_KERNEL::Exception("ParaDomainSelector::recvMesh : incoherent call in non_MPI mode");
591 // First stage : exchanging sizes
592 // ------------------------------
593 std::vector<int> tinyInfoDistant;
594 std::vector<std::string> tinyInfoLocalS;
595 std::vector<double> tinyInfoDistantD(1);
596 //Getting tiny info of local mesh to allow the distant proc to initialize and allocate
597 //the transmitted mesh.
600 MPI_Recv(&tinyVecSize, 1, MPI_INT, source, 1113, MPI_COMM_WORLD, &status);
601 tinyInfoDistant.resize(tinyVecSize);
602 std::fill(tinyInfoDistant.begin(),tinyInfoDistant.end(),0);
604 MPI_Recv(&tinyInfoDistant[0], tinyVecSize, MPI_INT,source,1112,MPI_COMM_WORLD, &status);
605 //there was tinyInfoLocal.push_back(mesh.getNumberOfCells());
606 int NumberOfCells=tinyInfoDistant[tinyVecSize-1];
609 MEDCoupling::DataArrayInt *v1Distant=MEDCoupling::DataArrayInt::New();
610 MEDCoupling::DataArrayDouble *v2Distant=MEDCoupling::DataArrayDouble::New();
611 //Building the right instance of copy of distant mesh.
612 MEDCoupling::MEDCouplingPointSet *distant_mesh_tmp=
613 MEDCoupling::MEDCouplingPointSet::BuildInstanceFromMeshType(
614 (MEDCoupling::MEDCouplingMeshType) tinyInfoDistant[0]);
615 std::vector<std::string> unusedTinyDistantSts;
616 mesh=dynamic_cast<MEDCoupling::MEDCouplingUMesh*> (distant_mesh_tmp);
618 mesh->resizeForUnserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
623 nbDistElem=v1Distant->getNbOfElems();
624 ptDist=v1Distant->getPointer();
626 MPI_Recv(ptDist, nbDistElem, MPI_INT, source,1111, MPI_COMM_WORLD, &status);
631 nbDistElem=v2Distant->getNbOfElems();
632 ptDist2=v2Distant->getPointer();
634 MPI_Recv(ptDist2, nbDistElem, MPI_DOUBLE,source, 1110, MPI_COMM_WORLD, &status);
635 //finish unserialization
636 mesh->unserialization(tinyInfoDistantD,tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
637 if(v1Distant) v1Distant->decrRef();
638 if(v2Distant) v2Distant->decrRef();
642 mesh=CreateEmptyMEDCouplingUMesh();
644 if (MyGlobals::_Verbose>600)
645 std::cout << "proc " << _rank << " : recvMesh '" << mesh->getName() << "' size " << mesh->getNumberOfCells() << " from " << source << std::endl;
649 #if !defined WIN32 && !defined __APPLE__
650 #include <sys/sysinfo.h>
654 * \brief Evaluate current memory usage and return the maximal one in KB
656 int MEDPARTITIONER::ParaDomainSelector::evaluateMemory() const
658 if ( _mesure_memory )
661 #if !defined WIN32 && !defined __APPLE__
663 int err = sysinfo( &si );
665 used_memory = (( si.totalram - si.freeram + si.totalswap - si.freeswap ) * si.mem_unit ) / 1024;
667 if ( used_memory > _max_memory )
668 _max_memory = used_memory;
671 _init_memory = used_memory;
673 return _max_memory - _init_memory;