1 // Copyright (C) 2007-2012 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.
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"
34 * \brief Constructor. Find out my rank and world size
36 MEDPARTITIONER::ParaDomainSelector::ParaDomainSelector(bool mesure_memory)
37 :_rank(0),_world_size(1), _nb_result_domains(-1), _init_time(0.0),
38 _mesure_memory(mesure_memory), _init_memory(0), _max_memory(0)
41 if (MyGlobals::_Rank==-1)
43 MPI_Init(0,0); //do once only
44 MPI_Comm_size(MPI_COMM_WORLD,&_world_size) ;
45 MPI_Comm_rank(MPI_COMM_WORLD,&_rank) ;
49 _world_size=MyGlobals::_World_Size;
50 _rank=MyGlobals::_Rank;
52 _init_time = MPI_Wtime();
57 if (MyGlobals::_Verbose>10)
58 std::cout << "WARNING : ParaDomainSelector contructor without parallel_mode World_Size=1 by default" << std::endl;
60 MyGlobals::_World_Size=_world_size;
61 MyGlobals::_Rank=_rank;
63 if (MyGlobals::_Verbose>200) std::cout << "proc " << MyGlobals::_Rank << " of " << MyGlobals::_World_Size << std::endl;
67 MEDPARTITIONER::ParaDomainSelector::~ParaDomainSelector()
72 * \brief Return true if is running on different hosts
74 bool MEDPARTITIONER::ParaDomainSelector::isOnDifferentHosts() const
77 if ( _world_size < 2 ) return false;
80 char name_here[ MPI_MAX_PROCESSOR_NAME+1 ], name_there[ MPI_MAX_PROCESSOR_NAME+1 ];
82 MPI_Get_processor_name( name_here, &size);
84 int next_proc = (rank() + 1) % nbProcs();
85 int prev_proc = (rank() - 1 + nbProcs()) % nbProcs();
89 MPI_Sendrecv((void*)&name_here[0], MPI_MAX_PROCESSOR_NAME, MPI_CHAR, next_proc, tag,
90 (void*)&name_there[0], MPI_MAX_PROCESSOR_NAME, MPI_CHAR, prev_proc, tag,
91 MPI_COMM_WORLD, &status);
93 //bug: (isOnDifferentHosts here and there) is not (isOnDifferentHosts somewhere)
94 //return string(name_here) != string(name_there);
98 if (std::string(name_here) != std::string(name_there))
100 MPI_Allreduce( &same, &sum_same, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );
101 return (sum_same != nbProcs());
107 * \brief Return true if the domain with domainIndex is to be loaded on this proc
108 * \param domainIndex - index of mesh domain
109 * \retval bool - to load or not
111 bool MEDPARTITIONER::ParaDomainSelector::isMyDomain(int domainIndex) const
114 return (_rank == getProcessorID( domainIndex ));
118 * \brief Return processor id where the domain with domainIndex resides
119 * \param domainIndex - index of mesh domain
120 * \retval int - processor id
122 int MEDPARTITIONER::ParaDomainSelector::getProcessorID(int domainIndex) const
125 return ( domainIndex % _world_size );
129 * \brief Gather info on nb of cell entities on each processor and return total nb.
132 * 1) for MED_CELL to know global id shift for domains at graph construction;
133 * 2) for sub-entity to know total nb of sub-entities before creating those of joints
135 void MEDPARTITIONER::ParaDomainSelector::gatherNbOf(const std::vector<ParaMEDMEM::MEDCouplingUMesh*>& domain_meshes)
138 // get nb of elems of each domain mesh
139 int nb_domains=domain_meshes.size();
140 std::vector<int> nb_elems(nb_domains*2, 0); //NumberOfCells & NumberOfNodes
141 for (int i=0; i<nb_domains; ++i)
142 if ( domain_meshes[i] )
144 nb_elems[i*2] = domain_meshes[i]->getNumberOfCells();
145 nb_elems[i*2+1] = domain_meshes[i]->getNumberOfNodes();
147 // receive nb of elems from other procs
148 std::vector<int> all_nb_elems;
149 if (MyGlobals::_World_Size==1)
151 all_nb_elems=nb_elems;
156 all_nb_elems.resize( nb_domains*2 );
157 MPI_Allreduce((void*)&nb_elems[0], (void*)&all_nb_elems[0], nb_domains*2, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
159 throw INTERP_KERNEL::Exception("not(HAVE_MPI2) incompatible with MPI_World_Size>1");
162 int total_nb_cells=0, total_nb_nodes=0;
163 for (int i=0; i<nb_domains; ++i)
165 total_nb_cells+=all_nb_elems[i*2];
166 total_nb_nodes+=all_nb_elems[i*2+1];
169 if (MyGlobals::_Is0verbose>10)
170 std::cout << "totalNbCells " << total_nb_cells << " totalNbNodes " << total_nb_nodes << std::endl;
172 std::vector<int>& cell_shift_by_domain=_cell_shift_by_domain;
173 std::vector<int>& node_shift_by_domain=_node_shift_by_domain;
174 std::vector<int>& face_shift_by_domain=_face_shift_by_domain;
176 std::vector< int > ordered_nbs_cell, ordered_nbs_node, domain_order(nb_domains);
177 ordered_nbs_cell.push_back(0);
178 ordered_nbs_node.push_back(0);
179 for (int iproc=0; iproc<nbProcs(); ++iproc)
180 for (int idomain=0; idomain<nb_domains; ++idomain)
181 if (getProcessorID( idomain )==iproc)
183 domain_order[idomain] = ordered_nbs_cell.size() - 1;
184 ordered_nbs_cell.push_back( ordered_nbs_cell.back() + all_nb_elems[idomain*2] );
185 ordered_nbs_node.push_back( ordered_nbs_node.back() + all_nb_elems[idomain*2+1] );
187 cell_shift_by_domain.resize( nb_domains+1, 0 );
188 node_shift_by_domain.resize( nb_domains+1, 0 );
189 face_shift_by_domain.resize( nb_domains+1, 0 );
190 for (int idomain=0; idomain<nb_domains; ++idomain)
192 cell_shift_by_domain[ idomain ] = ordered_nbs_cell[ domain_order[ idomain ]];
193 node_shift_by_domain[ idomain ] = ordered_nbs_node[ domain_order[ idomain ]];
195 cell_shift_by_domain.back() = ordered_nbs_cell.back(); // to know total nb of elements
196 node_shift_by_domain.back() = ordered_nbs_node.back(); // to know total nb of elements
198 if (MyGlobals::_Is0verbose>300)
200 std::cout << "proc " << MyGlobals::_Rank << " : cellShiftByDomain ";
201 for (int i=0; i<=nb_domains; ++i)
202 std::cout << cell_shift_by_domain[i] << "|";
203 std::cout << std::endl;
204 std::cout << "proc " << MyGlobals::_Rank << " : nodeShiftBy_domain ";
205 for (int i=0; i<=nb_domains; ++i)
206 std::cout << node_shift_by_domain[i] << "|";
207 std::cout << std::endl;
209 // fill _nb_vert_of_procs (is Vtxdist)
210 _nb_vert_of_procs.resize(_world_size+1);
211 _nb_vert_of_procs[0] = 0; // base = 0
212 for (int i=0; i<nb_domains; ++i)
214 int rankk = getProcessorID(i);
215 _nb_vert_of_procs[rankk+1] += all_nb_elems[i*2];
217 for (std::size_t i=1; i<_nb_vert_of_procs.size(); ++i)
218 _nb_vert_of_procs[i] += _nb_vert_of_procs[i-1]; // to CSR format : cumulated
220 if (MyGlobals::_Is0verbose>200)
222 std::cout << "proc " << MyGlobals::_Rank << " : gatherNbOf : vtxdist is ";
223 for (int i = 0; i <= _world_size; ++i)
224 std::cout << _nb_vert_of_procs[i] << " ";
225 std::cout << std::endl;
233 * \brief Return distribution of the graph vertices among the processors
234 * \retval int* - array containing nb of vertices (=cells) on all processors
236 * gatherNbOf() must be called before.
237 * The result array is to be used as the first arg of ParMETIS_V3_PartKway() and
238 * is freed by ParaDomainSelector.
240 int *MEDPARTITIONER::ParaDomainSelector::getProcVtxdist() const
243 if (_nb_vert_of_procs.empty())
244 throw INTERP_KERNEL::Exception("_nb_vert_of_procs not set");
245 return const_cast<int*>(& _nb_vert_of_procs[0]);
249 * \brief Return nb of cells in domains with lower index.
251 * gatherNbOf() must be called before.
252 * Result added to local id on given domain gives id in the whole distributed mesh
254 int MEDPARTITIONER::ParaDomainSelector::getDomainCellShift(int domainIndex) const
257 if (_cell_shift_by_domain.empty())
258 throw INTERP_KERNEL::Exception("_cell_shift_by_domain not set");
259 return _cell_shift_by_domain[domainIndex];
262 int MEDPARTITIONER::ParaDomainSelector::getDomainNodeShift(int domainIndex) const
265 if (_node_shift_by_domain.empty())
266 throw INTERP_KERNEL::Exception("_node_shift_by_domain not set");
267 return _node_shift_by_domain[domainIndex];
271 * \brief Return nb of nodes on processors with lower rank.
273 * gatherNbOf() must be called before.
274 * Result added to global id on this processor gives id in the whole distributed mesh
276 int MEDPARTITIONER::ParaDomainSelector::getProcNodeShift() const
279 if (_nb_vert_of_procs.empty())
280 throw INTERP_KERNEL::Exception("_nb_vert_of_procs not set");
281 return _nb_vert_of_procs[_rank];
285 * \brief Gather graphs from all processors into one
287 std::auto_ptr<MEDPARTITIONER::Graph> MEDPARTITIONER::ParaDomainSelector::gatherGraph(const Graph* graph) const
289 Graph* glob_graph = 0;
298 std::vector<int> index_size_of_proc( nbProcs() ); // index sizes - 1
299 for ( std::size_t i = 1; i < _nb_vert_of_procs.size(); ++i )
300 index_size_of_proc[i-1] = _nb_vert_of_procs[ i ] - _nb_vert_of_procs[ i-1 ];
302 int index_size = 1 + _cell_shift_by_domain.back();
303 int *graph_index = new int[ index_size ];
304 const int *index = graph->getGraph()->getIndex();
305 int *proc_index_displacement = const_cast<int*>( & _nb_vert_of_procs[0] );
307 MPI_Allgatherv((void*) (index+1), // send local index except first 0 (or 1)
308 index_size_of_proc[_rank], // index size on this proc
310 (void*) graph_index, // receive indices
311 & index_size_of_proc[0], // index size on each proc
312 proc_index_displacement, // displacement of each proc index
315 graph_index[0] = index[0]; // it is not overwritten thanks to proc_index_displacement[0]==1
317 // get sizes of graph values on each proc by the got indices of graphs
318 std::vector< int > value_size_of_proc( nbProcs() ), proc_value_displacement(1,0);
319 for ( int i = 0; i < nbProcs(); ++i )
321 if ( index_size_of_proc[i] > 0 )
322 value_size_of_proc[i] = graph_index[ proc_index_displacement[ i+1 ]-1 ] - graph_index[0];
324 value_size_of_proc[i] = 0;
325 proc_value_displacement.push_back( proc_value_displacement.back() + value_size_of_proc[i] );
328 // update graph_index
329 for ( int i = 1; i < nbProcs(); ++i )
331 int shift = graph_index[ proc_index_displacement[i]-1 ]-graph_index[0];
332 for ( int j = proc_index_displacement[i]; j < proc_index_displacement[i+1]; ++j )
333 graph_index[ j ] += shift;
340 int value_size = graph_index[ index_size-1 ] - graph_index[ 0 ];
341 int *graph_value = new int[ value_size ];
342 const int *value = graph->getGraph()->getValue();
344 MPI_Allgatherv((void*) value, // send local value
345 value_size_of_proc[_rank], // value size on this proc
347 (void*) graph_value, // receive values
348 & value_size_of_proc[0], // value size on each proc
349 & proc_value_displacement[0], // displacement of each proc value
357 int * partition = new int[ _cell_shift_by_domain.back() ];
358 const int* part = graph->getPart();
360 MPI_Allgatherv((void*) part, // send local partition
361 index_size_of_proc[_rank], // index size on this proc
363 (void*)(partition-1), // -1 compensates proc_index_displacement[0]==1
364 & index_size_of_proc[0], // index size on each proc
365 proc_index_displacement, // displacement of each proc index
373 // MEDPARTITIONER::SkyLineArray* array =
374 // new MEDPARTITIONER::SkyLineArray( index_size-1, value_size, graph_index, graph_value, true );
376 // glob_graph = new UserGraph( array, partition, index_size-1 );
384 return std::auto_ptr<Graph>( glob_graph );
389 * \brief Set nb of cell/cell pairs in a joint between domains
391 void MEDPARTITIONER::ParaDomainSelector::setNbCellPairs( int nb_cell_pairs, int dist_domain, int loc_domain )
393 // This method is needed for further computing global numbers of faces in joint.
394 // Store if both domains are on this proc else on one of procs only
395 if ( isMyDomain( dist_domain ) || dist_domain < loc_domain )
397 if ( _nb_cell_pairs_by_joint.empty() )
398 _nb_cell_pairs_by_joint.resize( _nb_result_domains*(_nb_result_domains+1), 0);
400 int joint_id = jointId( loc_domain, dist_domain );
401 _nb_cell_pairs_by_joint[ joint_id ] = nb_cell_pairs;
406 //================================================================================
408 * \brief Return nb of cell/cell pairs in a joint between domains on different procs
410 //================================================================================
412 int MEDPARTITIONER::ParaDomainSelector::getNbCellPairs( int dist_domain, int loc_domain ) const
416 int joint_id = jointId( loc_domain, dist_domain );
417 return _nb_cell_pairs_by_joint[ joint_id ];
420 //================================================================================
422 * \brief Gather size of each joint
424 //================================================================================
426 void MEDPARTITIONER::ParaDomainSelector::gatherNbCellPairs()
428 if ( _nb_cell_pairs_by_joint.empty() )
429 _nb_cell_pairs_by_joint.resize( _nb_result_domains*(_nb_result_domains+1), 0);
432 std::vector< int > send_buf = _nb_cell_pairs_by_joint;
434 MPI_Allreduce((void*)&send_buf[0],
435 (void*)&_nb_cell_pairs_by_joint[0],
436 _nb_cell_pairs_by_joint.size(),
437 MPI_INT, MPI_SUM, MPI_COMM_WORLD);
439 // check that the set nbs of cell pairs are correct,
440 // namely that each joint is treated on one proc only
441 for ( std::size_t j = 0; j < _nb_cell_pairs_by_joint.size(); ++j )
442 if ( _nb_cell_pairs_by_joint[j] != send_buf[j] && send_buf[j]>0 )
443 throw INTERP_KERNEL::Exception("invalid nb of cell pairs");
446 //================================================================================
448 * \brief Return the first global id of sub-entity for the joint
450 //================================================================================
452 int MEDPARTITIONER::ParaDomainSelector::getFisrtGlobalIdOfSubentity( int loc_domain, int dist_domain ) const
454 // total_nb_faces includes faces existing before creation of joint faces
455 // (got in gatherNbOf( MED_FACE )).
458 int total_nb_faces = _face_shift_by_domain.empty() ? 0 : _face_shift_by_domain.back();
459 int id = total_nb_faces + 1;
461 if ( _nb_cell_pairs_by_joint.empty() )
462 throw INTERP_KERNEL::Exception("gatherNbCellPairs() must be called before");
463 int joint_id = jointId( loc_domain, dist_domain );
464 for ( int j = 0; j < joint_id; ++j )
465 id += _nb_cell_pairs_by_joint[ j ];
470 //================================================================================
472 * \brief Send-receive local ids of joint faces
474 //================================================================================
476 int *MEDPARTITIONER::ParaDomainSelector::exchangeSubentityIds( int loc_domain, int dist_domain,
477 const std::vector<int>& loc_ids_here ) const
479 int* loc_ids_dist = new int[ loc_ids_here.size()];
481 int dest = getProcessorID( dist_domain );
482 int tag = 2002 + jointId( loc_domain, dist_domain );
484 MPI_Sendrecv((void*)&loc_ids_here[0], loc_ids_here.size(), MPI_INT, dest, tag,
485 (void*) loc_ids_dist, loc_ids_here.size(), MPI_INT, dest, tag,
486 MPI_COMM_WORLD, &status);
493 //================================================================================
495 * \brief Return identifier for a joint
497 //================================================================================
499 int MEDPARTITIONER::ParaDomainSelector::jointId( int local_domain, int distant_domain ) const
502 if (_nb_result_domains < 0)
503 throw INTERP_KERNEL::Exception("setNbDomains() must be called before");
505 if ( local_domain < distant_domain )
506 std::swap( local_domain, distant_domain );
507 return local_domain * _nb_result_domains + distant_domain;
511 //================================================================================
513 * \brief Return time passed from construction in seconds
515 //================================================================================
517 double MEDPARTITIONER::ParaDomainSelector::getPassedTime() const
520 return MPI_Wtime() - _init_time;
527 Sends content of \a mesh to processor \a target. To be used with \a recvMesh method.
528 \param mesh mesh to be sent
529 \param target processor id of the target
532 void MEDPARTITIONER::ParaDomainSelector::sendMesh(const ParaMEDMEM::MEDCouplingUMesh& mesh, int target) const
535 throw INTERP_KERNEL::Exception("ParaDomainSelector::sendMesh : incoherent call in non_MPI mode");
537 if (MyGlobals::_Verbose>600)
538 std::cout << "proc " << _rank << " : sendMesh '" << mesh.getName() << "' size " << mesh.getNumberOfCells() << " to " << target << std::endl;
539 // First stage : sending sizes
540 // ------------------------------
541 std::vector<int> tinyInfoLocal;
542 std::vector<std::string> tinyInfoLocalS;
543 std::vector<double> tinyInfoLocalD;
544 //Getting tiny info of local mesh to allow the distant proc to initialize and allocate
545 //the transmitted mesh.
546 mesh.getTinySerializationInformation(tinyInfoLocalD,tinyInfoLocal,tinyInfoLocalS);
547 tinyInfoLocal.push_back(mesh.getNumberOfCells());
548 int tinySize=tinyInfoLocal.size();
549 MPI_Send(&tinySize, 1, MPI_INT, target, 1113, MPI_COMM_WORLD);
550 MPI_Send(&tinyInfoLocal[0], tinyInfoLocal.size(), MPI_INT, target, 1112, MPI_COMM_WORLD);
552 if (mesh.getNumberOfCells()>0) //no sends if empty
554 ParaMEDMEM::DataArrayInt *v1Local=0;
555 ParaMEDMEM::DataArrayDouble *v2Local=0;
556 //serialization of local mesh to send data to distant proc.
557 mesh.serialize(v1Local,v2Local);
560 if(v1Local) //if empty getNbOfElems() is 1!
562 nbLocalElems=v1Local->getNbOfElems(); // if empty be 1!
563 ptLocal=v1Local->getPointer();
565 MPI_Send(ptLocal, nbLocalElems, MPI_INT, target, 1111, MPI_COMM_WORLD);
568 if(v2Local) //if empty be 0!
570 nbLocalElems2=v2Local->getNbOfElems();
571 ptLocal2=v2Local->getPointer();
573 MPI_Send(ptLocal2, nbLocalElems2, MPI_DOUBLE, target, 1110, MPI_COMM_WORLD);
574 if(v1Local) v1Local->decrRef();
575 if(v2Local) v2Local->decrRef();
580 /*! Receives messages from proc \a source to fill mesh \a mesh.
581 To be used with \a sendMesh method.
582 \param mesh pointer to mesh that is filled
583 \param source processor id of the incoming messages
585 void MEDPARTITIONER::ParaDomainSelector::recvMesh(ParaMEDMEM::MEDCouplingUMesh*& mesh, int source)const
588 throw INTERP_KERNEL::Exception("ParaDomainSelector::recvMesh : incoherent call in non_MPI mode");
590 // First stage : exchanging sizes
591 // ------------------------------
592 std::vector<int> tinyInfoDistant;
593 std::vector<std::string> tinyInfoLocalS;
594 std::vector<double> tinyInfoDistantD(1);
595 //Getting tiny info of local mesh to allow the distant proc to initialize and allocate
596 //the transmitted mesh.
599 MPI_Recv(&tinyVecSize, 1, MPI_INT, source, 1113, MPI_COMM_WORLD, &status);
600 tinyInfoDistant.resize(tinyVecSize);
601 std::fill(tinyInfoDistant.begin(),tinyInfoDistant.end(),0);
603 MPI_Recv(&tinyInfoDistant[0], tinyVecSize, MPI_INT,source,1112,MPI_COMM_WORLD, &status);
604 //there was tinyInfoLocal.push_back(mesh.getNumberOfCells());
605 int NumberOfCells=tinyInfoDistant[tinyVecSize-1];
608 ParaMEDMEM::DataArrayInt *v1Distant=ParaMEDMEM::DataArrayInt::New();
609 ParaMEDMEM::DataArrayDouble *v2Distant=ParaMEDMEM::DataArrayDouble::New();
610 //Building the right instance of copy of distant mesh.
611 ParaMEDMEM::MEDCouplingPointSet *distant_mesh_tmp=
612 ParaMEDMEM::MEDCouplingPointSet::BuildInstanceFromMeshType(
613 (ParaMEDMEM::MEDCouplingMeshType) tinyInfoDistant[0]);
614 std::vector<std::string> unusedTinyDistantSts;
615 mesh=dynamic_cast<ParaMEDMEM::MEDCouplingUMesh*> (distant_mesh_tmp);
617 mesh->resizeForUnserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
622 nbDistElem=v1Distant->getNbOfElems();
623 ptDist=v1Distant->getPointer();
625 MPI_Recv(ptDist, nbDistElem, MPI_INT, source,1111, MPI_COMM_WORLD, &status);
630 nbDistElem=v2Distant->getNbOfElems();
631 ptDist2=v2Distant->getPointer();
633 MPI_Recv(ptDist2, nbDistElem, MPI_DOUBLE,source, 1110, MPI_COMM_WORLD, &status);
634 //finish unserialization
635 mesh->unserialization(tinyInfoDistantD,tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
636 if(v1Distant) v1Distant->decrRef();
637 if(v2Distant) v2Distant->decrRef();
641 mesh=CreateEmptyMEDCouplingUMesh();
643 if (MyGlobals::_Verbose>600)
644 std::cout << "proc " << _rank << " : recvMesh '" << mesh->getName() << "' size " << mesh->getNumberOfCells() << " from " << source << std::endl;
649 #include <sys/sysinfo.h>
653 * \brief Evaluate current memory usage and return the maximal one in KB
655 int MEDPARTITIONER::ParaDomainSelector::evaluateMemory() const
657 if ( _mesure_memory )
662 int err = sysinfo( &si );
664 used_memory = (( si.totalram - si.freeram + si.totalswap - si.freeswap ) * si.mem_unit ) / 1024;
666 if ( used_memory > _max_memory )
667 _max_memory = used_memory;
670 _init_memory = used_memory;
672 return _max_memory - _init_memory;