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());
106 * \brief Return true if the domain with domainIndex is to be loaded on this proc
107 * \param domainIndex - index of mesh domain
108 * \retval bool - to load or not
110 bool MEDPARTITIONER::ParaDomainSelector::isMyDomain(int domainIndex) const
113 return (_rank == getProcessorID( domainIndex ));
117 * \brief Return processor id where the domain with domainIndex resides
118 * \param domainIndex - index of mesh domain
119 * \retval int - processor id
121 int MEDPARTITIONER::ParaDomainSelector::getProcessorID(int domainIndex) const
124 return ( domainIndex % _world_size );
128 * \brief Gather info on nb of cell entities on each processor and return total nb.
131 * 1) for MED_CELL to know global id shift for domains at graph construction;
132 * 2) for sub-entity to know total nb of sub-entities before creating those of joints
134 void MEDPARTITIONER::ParaDomainSelector::gatherNbOf(const std::vector<ParaMEDMEM::MEDCouplingUMesh*>& domain_meshes)
137 // get nb of elems of each domain mesh
138 int nb_domains=domain_meshes.size();
139 std::vector<int> nb_elems(nb_domains*2, 0); //NumberOfCells & NumberOfNodes
140 for (int i=0; i<nb_domains; ++i)
141 if ( domain_meshes[i] )
143 nb_elems[i*2] = domain_meshes[i]->getNumberOfCells();
144 nb_elems[i*2+1] = domain_meshes[i]->getNumberOfNodes();
146 // receive nb of elems from other procs
147 std::vector<int> all_nb_elems;
148 if (MyGlobals::_World_Size==1)
150 all_nb_elems=nb_elems;
155 all_nb_elems.resize( nb_domains*2 );
156 MPI_Allreduce((void*)&nb_elems[0], (void*)&all_nb_elems[0], nb_domains*2, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
158 throw INTERP_KERNEL::Exception("not(HAVE_MPI2) incompatible with MPI_World_Size>1");
161 int total_nb_cells=0, total_nb_nodes=0;
162 for (int i=0; i<nb_domains; ++i)
164 total_nb_cells+=all_nb_elems[i*2];
165 total_nb_nodes+=all_nb_elems[i*2+1];
168 if (MyGlobals::_Is0verbose>10)
169 std::cout << "totalNbCells " << total_nb_cells << " totalNbNodes " << total_nb_nodes << std::endl;
171 std::vector<int>& cell_shift_by_domain=_cell_shift_by_domain;
172 std::vector<int>& node_shift_by_domain=_node_shift_by_domain;
173 std::vector<int>& face_shift_by_domain=_face_shift_by_domain;
175 std::vector< int > ordered_nbs_cell, ordered_nbs_node, domain_order(nb_domains);
176 ordered_nbs_cell.push_back(0);
177 ordered_nbs_node.push_back(0);
178 for (int iproc=0; iproc<nbProcs(); ++iproc)
179 for (int idomain=0; idomain<nb_domains; ++idomain)
180 if (getProcessorID( idomain )==iproc)
182 domain_order[idomain] = ordered_nbs_cell.size() - 1;
183 ordered_nbs_cell.push_back( ordered_nbs_cell.back() + all_nb_elems[idomain*2] );
184 ordered_nbs_node.push_back( ordered_nbs_node.back() + all_nb_elems[idomain*2+1] );
186 cell_shift_by_domain.resize( nb_domains+1, 0 );
187 node_shift_by_domain.resize( nb_domains+1, 0 );
188 face_shift_by_domain.resize( nb_domains+1, 0 );
189 for (int idomain=0; idomain<nb_domains; ++idomain)
191 cell_shift_by_domain[ idomain ] = ordered_nbs_cell[ domain_order[ idomain ]];
192 node_shift_by_domain[ idomain ] = ordered_nbs_node[ domain_order[ idomain ]];
194 cell_shift_by_domain.back() = ordered_nbs_cell.back(); // to know total nb of elements
195 node_shift_by_domain.back() = ordered_nbs_node.back(); // to know total nb of elements
197 if (MyGlobals::_Is0verbose>300)
199 std::cout << "proc " << MyGlobals::_Rank << " : cellShiftByDomain ";
200 for (int i=0; i<=nb_domains; ++i)
201 std::cout << cell_shift_by_domain[i] << "|";
202 std::cout << std::endl;
203 std::cout << "proc " << MyGlobals::_Rank << " : nodeShiftBy_domain ";
204 for (int i=0; i<=nb_domains; ++i)
205 std::cout << node_shift_by_domain[i] << "|";
206 std::cout << std::endl;
208 // fill _nb_vert_of_procs (is Vtxdist)
209 _nb_vert_of_procs.resize(_world_size+1);
210 _nb_vert_of_procs[0] = 0; // base = 0
211 for (int i=0; i<nb_domains; ++i)
213 int rankk = getProcessorID(i);
214 _nb_vert_of_procs[rankk+1] += all_nb_elems[i*2];
216 for (std::size_t i=1; i<_nb_vert_of_procs.size(); ++i)
217 _nb_vert_of_procs[i] += _nb_vert_of_procs[i-1]; // to CSR format : cumulated
219 if (MyGlobals::_Is0verbose>200)
221 std::cout << "proc " << MyGlobals::_Rank << " : gatherNbOf : vtxdist is ";
222 for (int i = 0; i <= _world_size; ++i)
223 std::cout << _nb_vert_of_procs[i] << " ";
224 std::cout << std::endl;
232 * \brief Return distribution of the graph vertices among the processors
233 * \retval int* - array containing nb of vertices (=cells) on all processors
235 * gatherNbOf() must be called before.
236 * The result array is to be used as the first arg of ParMETIS_V3_PartKway() and
237 * is freed by ParaDomainSelector.
239 int *MEDPARTITIONER::ParaDomainSelector::getProcVtxdist() const
242 if (_nb_vert_of_procs.empty())
243 throw INTERP_KERNEL::Exception("_nb_vert_of_procs not set");
244 return const_cast<int*>(& _nb_vert_of_procs[0]);
248 * \brief Return nb of cells in domains with lower index.
250 * gatherNbOf() must be called before.
251 * Result added to local id on given domain gives id in the whole distributed mesh
253 int MEDPARTITIONER::ParaDomainSelector::getDomainCellShift(int domainIndex) const
256 if (_cell_shift_by_domain.empty())
257 throw INTERP_KERNEL::Exception("_cell_shift_by_domain not set");
258 return _cell_shift_by_domain[domainIndex];
261 int MEDPARTITIONER::ParaDomainSelector::getDomainNodeShift(int domainIndex) const
264 if (_node_shift_by_domain.empty())
265 throw INTERP_KERNEL::Exception("_node_shift_by_domain not set");
266 return _node_shift_by_domain[domainIndex];
270 * \brief Return nb of nodes on processors with lower rank.
272 * gatherNbOf() must be called before.
273 * Result added to global id on this processor gives id in the whole distributed mesh
275 int MEDPARTITIONER::ParaDomainSelector::getProcNodeShift() const
278 if (_nb_vert_of_procs.empty())
279 throw INTERP_KERNEL::Exception("_nb_vert_of_procs not set");
280 return _nb_vert_of_procs[_rank];
284 * \brief Gather graphs from all processors into one
286 std::auto_ptr<MEDPARTITIONER::Graph> MEDPARTITIONER::ParaDomainSelector::gatherGraph(const Graph* graph) const
288 Graph* glob_graph = 0;
297 std::vector<int> index_size_of_proc( nbProcs() ); // index sizes - 1
298 for ( std::size_t i = 1; i < _nb_vert_of_procs.size(); ++i )
299 index_size_of_proc[i-1] = _nb_vert_of_procs[ i ] - _nb_vert_of_procs[ i-1 ];
301 int index_size = 1 + _cell_shift_by_domain.back();
302 int *graph_index = new int[ index_size ];
303 const int *index = graph->getGraph()->getIndex();
304 int *proc_index_displacement = const_cast<int*>( & _nb_vert_of_procs[0] );
306 MPI_Allgatherv((void*) (index+1), // send local index except first 0 (or 1)
307 index_size_of_proc[_rank], // index size on this proc
309 (void*) graph_index, // receive indices
310 & index_size_of_proc[0], // index size on each proc
311 proc_index_displacement, // displacement of each proc index
314 graph_index[0] = index[0]; // it is not overwritten thanks to proc_index_displacement[0]==1
316 // get sizes of graph values on each proc by the got indices of graphs
317 std::vector< int > value_size_of_proc( nbProcs() ), proc_value_displacement(1,0);
318 for ( int i = 0; i < nbProcs(); ++i )
320 if ( index_size_of_proc[i] > 0 )
321 value_size_of_proc[i] = graph_index[ proc_index_displacement[ i+1 ]-1 ] - graph_index[0];
323 value_size_of_proc[i] = 0;
324 proc_value_displacement.push_back( proc_value_displacement.back() + value_size_of_proc[i] );
327 // update graph_index
328 for ( int i = 1; i < nbProcs(); ++i )
330 int shift = graph_index[ proc_index_displacement[i]-1 ]-graph_index[0];
331 for ( int j = proc_index_displacement[i]; j < proc_index_displacement[i+1]; ++j )
332 graph_index[ j ] += shift;
339 int value_size = graph_index[ index_size-1 ] - graph_index[ 0 ];
340 int *graph_value = new int[ value_size ];
341 const int *value = graph->getGraph()->getValue();
343 MPI_Allgatherv((void*) value, // send local value
344 value_size_of_proc[_rank], // value size on this proc
346 (void*) graph_value, // receive values
347 & value_size_of_proc[0], // value size on each proc
348 & proc_value_displacement[0], // displacement of each proc value
356 int * partition = new int[ _cell_shift_by_domain.back() ];
357 const int* part = graph->getPart();
359 MPI_Allgatherv((void*) part, // send local partition
360 index_size_of_proc[_rank], // index size on this proc
362 (void*)(partition-1), // -1 compensates proc_index_displacement[0]==1
363 & index_size_of_proc[0], // index size on each proc
364 proc_index_displacement, // displacement of each proc index
372 // MEDPARTITIONER::SkyLineArray* array =
373 // new MEDPARTITIONER::SkyLineArray( index_size-1, value_size, graph_index, graph_value, true );
375 // glob_graph = new UserGraph( array, partition, index_size-1 );
383 return std::auto_ptr<Graph>( glob_graph );
388 * \brief Set nb of cell/cell pairs in a joint between domains
390 void MEDPARTITIONER::ParaDomainSelector::setNbCellPairs( int nb_cell_pairs, int dist_domain, int loc_domain )
392 // This method is needed for further computing global numbers of faces in joint.
393 // Store if both domains are on this proc else on one of procs only
394 if ( isMyDomain( dist_domain ) || dist_domain < loc_domain )
396 if ( _nb_cell_pairs_by_joint.empty() )
397 _nb_cell_pairs_by_joint.resize( _nb_result_domains*(_nb_result_domains+1), 0);
399 int joint_id = jointId( loc_domain, dist_domain );
400 _nb_cell_pairs_by_joint[ joint_id ] = nb_cell_pairs;
405 //================================================================================
407 * \brief Return nb of cell/cell pairs in a joint between domains on different procs
409 //================================================================================
411 int MEDPARTITIONER::ParaDomainSelector::getNbCellPairs( int dist_domain, int loc_domain ) const
415 int joint_id = jointId( loc_domain, dist_domain );
416 return _nb_cell_pairs_by_joint[ joint_id ];
419 //================================================================================
421 * \brief Gather size of each joint
423 //================================================================================
425 void MEDPARTITIONER::ParaDomainSelector::gatherNbCellPairs()
427 if ( _nb_cell_pairs_by_joint.empty() )
428 _nb_cell_pairs_by_joint.resize( _nb_result_domains*(_nb_result_domains+1), 0);
431 std::vector< int > send_buf = _nb_cell_pairs_by_joint;
433 MPI_Allreduce((void*)&send_buf[0],
434 (void*)&_nb_cell_pairs_by_joint[0],
435 _nb_cell_pairs_by_joint.size(),
436 MPI_INT, MPI_SUM, MPI_COMM_WORLD);
438 // check that the set nbs of cell pairs are correct,
439 // namely that each joint is treated on one proc only
440 for ( std::size_t j = 0; j < _nb_cell_pairs_by_joint.size(); ++j )
441 if ( _nb_cell_pairs_by_joint[j] != send_buf[j] && send_buf[j]>0 )
442 throw INTERP_KERNEL::Exception("invalid nb of cell pairs");
445 //================================================================================
447 * \brief Return the first global id of sub-entity for the joint
449 //================================================================================
451 int MEDPARTITIONER::ParaDomainSelector::getFisrtGlobalIdOfSubentity( int loc_domain, int dist_domain ) const
453 // total_nb_faces includes faces existing before creation of joint faces
454 // (got in gatherNbOf( MED_FACE )).
457 int total_nb_faces = _face_shift_by_domain.empty() ? 0 : _face_shift_by_domain.back();
458 int id = total_nb_faces + 1;
460 if ( _nb_cell_pairs_by_joint.empty() )
461 throw INTERP_KERNEL::Exception("gatherNbCellPairs() must be called before");
462 int joint_id = jointId( loc_domain, dist_domain );
463 for ( int j = 0; j < joint_id; ++j )
464 id += _nb_cell_pairs_by_joint[ j ];
469 //================================================================================
471 * \brief Send-receive local ids of joint faces
473 //================================================================================
475 int *MEDPARTITIONER::ParaDomainSelector::exchangeSubentityIds( int loc_domain, int dist_domain,
476 const std::vector<int>& loc_ids_here ) const
478 int* loc_ids_dist = new int[ loc_ids_here.size()];
480 int dest = getProcessorID( dist_domain );
481 int tag = 2002 + jointId( loc_domain, dist_domain );
483 MPI_Sendrecv((void*)&loc_ids_here[0], loc_ids_here.size(), MPI_INT, dest, tag,
484 (void*) loc_ids_dist, loc_ids_here.size(), MPI_INT, dest, tag,
485 MPI_COMM_WORLD, &status);
492 //================================================================================
494 * \brief Return identifier for a joint
496 //================================================================================
498 int MEDPARTITIONER::ParaDomainSelector::jointId( int local_domain, int distant_domain ) const
501 if (_nb_result_domains < 0)
502 throw INTERP_KERNEL::Exception("setNbDomains() must be called before");
504 if ( local_domain < distant_domain )
505 std::swap( local_domain, distant_domain );
506 return local_domain * _nb_result_domains + distant_domain;
510 //================================================================================
512 * \brief Return time passed from construction in seconds
514 //================================================================================
516 double MEDPARTITIONER::ParaDomainSelector::getPassedTime() const
519 return MPI_Wtime() - _init_time;
526 Sends content of \a mesh to processor \a target. To be used with \a recvMesh method.
527 \param mesh mesh to be sent
528 \param target processor id of the target
531 void MEDPARTITIONER::ParaDomainSelector::sendMesh(const ParaMEDMEM::MEDCouplingUMesh& mesh, int target) const
534 throw INTERP_KERNEL::Exception("ParaDomainSelector::sendMesh : incoherent call in non_MPI mode");
536 if (MyGlobals::_Verbose>600)
537 std::cout << "proc " << _rank << " : sendMesh '" << mesh.getName() << "' size " << mesh.getNumberOfCells() << " to " << target << std::endl;
538 // First stage : sending sizes
539 // ------------------------------
540 std::vector<int> tinyInfoLocal;
541 std::vector<std::string> tinyInfoLocalS;
542 std::vector<double> tinyInfoLocalD;
543 //Getting tiny info of local mesh to allow the distant proc to initialize and allocate
544 //the transmitted mesh.
545 mesh.getTinySerializationInformation(tinyInfoLocalD,tinyInfoLocal,tinyInfoLocalS);
546 tinyInfoLocal.push_back(mesh.getNumberOfCells());
547 int tinySize=tinyInfoLocal.size();
548 MPI_Send(&tinySize, 1, MPI_INT, target, 1113, MPI_COMM_WORLD);
549 MPI_Send(&tinyInfoLocal[0], tinyInfoLocal.size(), MPI_INT, target, 1112, MPI_COMM_WORLD);
551 if (mesh.getNumberOfCells()>0) //no sends if empty
553 ParaMEDMEM::DataArrayInt *v1Local=0;
554 ParaMEDMEM::DataArrayDouble *v2Local=0;
555 //serialization of local mesh to send data to distant proc.
556 mesh.serialize(v1Local,v2Local);
559 if(v1Local) //if empty getNbOfElems() is 1!
561 nbLocalElems=v1Local->getNbOfElems(); // if empty be 1!
562 ptLocal=v1Local->getPointer();
564 MPI_Send(ptLocal, nbLocalElems, MPI_INT, target, 1111, MPI_COMM_WORLD);
567 if(v2Local) //if empty be 0!
569 nbLocalElems2=v2Local->getNbOfElems();
570 ptLocal2=v2Local->getPointer();
572 MPI_Send(ptLocal2, nbLocalElems2, MPI_DOUBLE, target, 1110, MPI_COMM_WORLD);
573 if(v1Local) v1Local->decrRef();
574 if(v2Local) v2Local->decrRef();
579 /*! Receives messages from proc \a source to fill mesh \a mesh.
580 To be used with \a sendMesh method.
581 \param mesh pointer to mesh that is filled
582 \param source processor id of the incoming messages
584 void MEDPARTITIONER::ParaDomainSelector::recvMesh(ParaMEDMEM::MEDCouplingUMesh*& mesh, int source)const
587 throw INTERP_KERNEL::Exception("ParaDomainSelector::recvMesh : incoherent call in non_MPI mode");
589 // First stage : exchanging sizes
590 // ------------------------------
591 std::vector<int> tinyInfoDistant;
592 std::vector<std::string> tinyInfoLocalS;
593 std::vector<double> tinyInfoDistantD(1);
594 //Getting tiny info of local mesh to allow the distant proc to initialize and allocate
595 //the transmitted mesh.
598 MPI_Recv(&tinyVecSize, 1, MPI_INT, source, 1113, MPI_COMM_WORLD, &status);
599 tinyInfoDistant.resize(tinyVecSize);
600 std::fill(tinyInfoDistant.begin(),tinyInfoDistant.end(),0);
602 MPI_Recv(&tinyInfoDistant[0], tinyVecSize, MPI_INT,source,1112,MPI_COMM_WORLD, &status);
603 //there was tinyInfoLocal.push_back(mesh.getNumberOfCells());
604 int NumberOfCells=tinyInfoDistant[tinyVecSize-1];
607 ParaMEDMEM::DataArrayInt *v1Distant=ParaMEDMEM::DataArrayInt::New();
608 ParaMEDMEM::DataArrayDouble *v2Distant=ParaMEDMEM::DataArrayDouble::New();
609 //Building the right instance of copy of distant mesh.
610 ParaMEDMEM::MEDCouplingPointSet *distant_mesh_tmp=
611 ParaMEDMEM::MEDCouplingPointSet::BuildInstanceFromMeshType(
612 (ParaMEDMEM::MEDCouplingMeshType) tinyInfoDistant[0]);
613 std::vector<std::string> unusedTinyDistantSts;
614 mesh=dynamic_cast<ParaMEDMEM::MEDCouplingUMesh*> (distant_mesh_tmp);
616 mesh->resizeForUnserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
621 nbDistElem=v1Distant->getNbOfElems();
622 ptDist=v1Distant->getPointer();
624 MPI_Recv(ptDist, nbDistElem, MPI_INT, source,1111, MPI_COMM_WORLD, &status);
629 nbDistElem=v2Distant->getNbOfElems();
630 ptDist2=v2Distant->getPointer();
632 MPI_Recv(ptDist2, nbDistElem, MPI_DOUBLE,source, 1110, MPI_COMM_WORLD, &status);
633 //finish unserialization
634 mesh->unserialization(tinyInfoDistantD,tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts);
635 if(v1Distant) v1Distant->decrRef();
636 if(v2Distant) v2Distant->decrRef();
640 mesh=CreateEmptyMEDCouplingUMesh();
642 if (MyGlobals::_Verbose>600)
643 std::cout << "proc " << _rank << " : recvMesh '" << mesh->getName() << "' size " << mesh->getNumberOfCells() << " from " << source << std::endl;
648 #include <sys/sysinfo.h>
652 * \brief Evaluate current memory usage and return the maximal one in KB
654 int MEDPARTITIONER::ParaDomainSelector::evaluateMemory() const
656 if ( _mesure_memory )
661 int err = sysinfo( &si );
663 used_memory = (( si.totalram - si.freeram + si.totalswap - si.freeswap ) * si.mem_unit ) / 1024;
665 if ( used_memory > _max_memory )
666 _max_memory = used_memory;
669 _init_memory = used_memory;
671 return _max_memory - _init_memory;