Mechanics of Defect Evolution Library
Loading...
Searching...
No Matches
ParticleSystemParallel.h
Go to the documentation of this file.
1/* This file is part of MoDELib, the Mechanics Of Defects Evolution Library.
2 *
3 *
4 * MODEL is distributed without any warranty under the
5 * GNU General Public License (GPL) v2 <http://www.gnu.org/licenses/>.
6 */
7
8#ifndef _MODEL_MPI_
9#error YOU NEED TO DEFINE _MODEL_MPI_ TO USE THIS HEADER
10#endif
11
12#ifndef _MODEL_ParticleSystemParallel_h_
13#define _MODEL_ParticleSystemParallel_h_
14
15#ifdef _MODEL_MPI_
16#include <ModelMPIbase.h>
17//#include <metis.h> // partitioner
18#endif
19
20#ifdef _OPENMP
21#include <omp.h>
22#endif
23
24
25#include <vector>
26#include <time.h>
27
28#include <modelMacros.h> // model_checkInput
29#include <SequentialOutputFile.h>
31#include <PointSource.h>
32#include <FieldPoint.h>
33#include <ParticleSystemBase.h>
34#include <LPTpartitioner.h>
35
36
37namespace model {
38
39 /**************************************************************************/
40 /**************************************************************************/
43 template <typename _ParticleType>
45 /* inheritance */ public ParticleSystemBase<_ParticleType>,
46 /* inheritance */ public ModelMPIbase
47 {
48
49 typedef ParticleSystemBase<_ParticleType> ParticleSystemBaseType;
50 typedef _ParticleType ParticleType; // make ParticleType available outside the class
51 typedef typename ParticleSystemBaseType::SpatialCellObserverType SpatialCellObserverType;
52
53 public:
54
57
58
59 /**********************************************************************/
61 /* init list */ ModelMPIbase(argc,argv)
62 {
63
64 }
65
66 /**********************************************************************/
68 {
69
70 }
71
72 /**********************************************************************/
73 template <typename ...AdditionalConstructorTypes>
74 ParticleType* addParticle(const typename ParticleSystemBaseType::PositionType& p, const AdditionalConstructorTypes&... args)
75 {
78 assert(mpiInitialized() && "CANNOT INSERT PARTICLES BEFORE CALLING ModelMPIbase::initMPI(argc,argv).");
79 // partitionIsValid=false;
81 }
82
83
84 // /**********************************************************************/
85 // template <typename T>
86 // void partionSystem(LPTpartitioner<T>& partitioner)
87 // {/*! @param[in] useCellPartitioner Flag that determines whether the
88 // * ParticleSystem is partitioned by particles or by cells.
89 // *
90 // * Partitions particles among processors.
91 // */
92 //
93 //
94 // // Clear and call the LPT partitioner
95 // partitioner.clear();
96 //
97 // for (typename ParticleContainerType::iterator pIter =this->begin();
98 // /* */ pIter!=this->end();
99 // /* */ ++pIter)
100 // {
101 // partitioner.insert(pIter->pCell->n2Weight(),&*pIter);
102 // }
103 //
104 // partitioner.partition(this->mpiProcs());
105 // }
106
107 /**********************************************************************/
108 template <typename FieldType,typename... OtherSourceTypes>
109 void computeNeighborField(const OtherSourceTypes&... otherSources)
110 {
115
116 MPI_Barrier(MPI_COMM_WORLD);
117
118
121 typedef LPTpartitioner<ParticleType> PartitionerType;
122 PartitionerType lpt;
123
125 for (typename ParticleContainerType::iterator pIter =this->begin();
126 /* */ pIter!=this->end();
127 /* */ ++pIter)
128 {
129 //lpt.insert(pIter->pCell->n2Weight(),&*pIter);
130 lpt.insert(pIter->pCell->neighborSize(),&*pIter);
131 }
133 lpt.partition(this->mpiProcs());
134
135 // Resize FieldPointBase<ParticleType,FieldType>::resultVector
137
138 // Assign mpiID
139 size_t binOffset(0);
140 for (size_t b=0;b<lpt.bins();++b)
141 {
142 const size_t binSize(lpt.bin(b).size());
143 for (size_t p=0;p<binSize;++p)
144 {
145 lpt.bin(b)[p]->set_mpiID(binOffset+p);
146 }
147 binOffset+=binSize;
148 }
149
150
151#ifdef _OPENMP
152#pragma omp parallel for
153#endif
154 for (unsigned int k=0; k<lpt.bin(this->mpiRank()).size();++k)
155 {
156 // Nearest-neighbor interaction
157 static_cast<FieldPointBase<ParticleType,FieldType>* const>(lpt.bin(this->mpiRank())[k]) -> setZero();
158
159 // loop over neighbor cells of the current particle
160 for (typename SpatialCellType::CellMapType::const_iterator cIter=lpt.bin(this->mpiRank())[k]->neighborCellsBegin();cIter!=lpt.bin(this->mpiRank())[k]->neighborCellsEnd();++cIter)
161 {
162 // loop over source particles in the neighbor cell
163 for(typename SpatialCellType::ParticleContainerType::const_iterator sIter=cIter->second->particleBegin();sIter!=cIter->second->particleEnd();++sIter) // loop over neighbor particles
164 {
165 // lpt.bin(this->mpiRank())[k]->template field<FieldType>() += FieldType::compute(**sIter,*lpt.bin(this->mpiRank())[k]);
166 if(static_cast<const SingleSourcePoint<ParticleType,FieldType>* const>(*sIter)->enabled)
167 {
168 *static_cast<FieldPointBase<ParticleType,FieldType>* const>(lpt.bin(this->mpiRank())[k]) += FieldType::compute(**sIter,*lpt.bin(this->mpiRank())[k]);
169 }
170 }
171 }
172
173 // Non nearest-neighbor interaction
174 if(FieldType::use_multipole)
175 {
176 typename SpatialCellType::CellMapType farCells(lpt.bin(this->mpiRank())[k]->template farCells<_ParticleType>());
177 // lpt.bin(this->mpiRank())[k]->template field<FieldType>() += FieldType::multipole(*lpt.bin(this->mpiRank())[k],farCells);
178 *static_cast<FieldPointBase<ParticleType,FieldType>* const>(lpt.bin(this->mpiRank())[k]) += FieldType::multipole(*lpt.bin(this->mpiRank())[k],farCells);
179 }
180
181 // Add contribution of other sources
182 if(sizeof...(OtherSourceTypes))
183 {
184 *static_cast<FieldPointBase<ParticleType,FieldType>* const>(lpt.bin(this->mpiRank())[k]) += FieldType::addSourceContribution(*lpt.bin(this->mpiRank())[k],otherSources...);
185 }
186 }
187
189 std::vector<int> interactionSizeVector(this->mpiProcs());
190 std::vector<int> interactionRankOffsetVector(this->mpiProcs());
191
192 size_t dofOffset(0);
193 for (int k=0;k<this->mpiProcs();++k)
194 {
195 // Store the number
196 interactionSizeVector[k]=lpt.bin(k).size()*FieldPointBase<ParticleType,FieldType>::DataPerParticle;
197 interactionRankOffsetVector[k]=dofOffset;
199 }
200
201 MPI_Allgatherv(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,
203 &interactionSizeVector[0],&interactionRankOffsetVector[0],
204 MPI_DOUBLE,MPI_COMM_WORLD);
205 }
206
207
208 /**********************************************************************/
209 template <typename OtherParticleType, typename FieldType,typename... OtherSourceTypes>
210 void computeField(std::deque<OtherParticleType>& fpDeq,const OtherSourceTypes&... otherSources) const
211 {
214
215 MPI_Barrier(MPI_COMM_WORLD);
216
217 // LPTpartitioner<ParticleType> lpt;
218 typedef LPTpartitioner<OtherParticleType> PartitionerType;
219 PartitionerType lpt;
220 // partionSystem(lpt);
221
222 for (typename std::deque<OtherParticleType>::iterator pIter =fpDeq.begin();
223 /* */ pIter!=fpDeq.end();
224 /* */ ++pIter)
225 {
226 typename SpatialCellType::CellMapType cellMap(pIter->template neighborCells<_ParticleType>());
227 double weight(0.0);
228 for (typename SpatialCellType::CellMapType::const_iterator cIter=cellMap.begin();cIter!=cellMap.end();++cIter)
229 {
230 weight += cIter->second->size();
231 }
232
233 lpt.insert(weight,&*pIter);
234 }
235
236 lpt.partition(this->mpiProcs());
237
238
239
240 // Resize FieldPointBase<OtherParticleType,FieldType>::resultVector
242
243 // Assign mpiID
244 size_t binOffset(0);
245 for (size_t b=0;b<lpt.bins();++b)
246 {
247 const size_t binSize(lpt.bin(b).size());
248 for (size_t p=0;p<binSize;++p)
249 {
250 lpt.bin(b)[p]->set_mpiID(binOffset+p);
251 }
252 binOffset+=binSize;
253 }
254
255
256#ifdef _OPENMP
257#pragma omp parallel for
258#endif
259 for (unsigned int k=0; k<lpt.bin(this->mpiRank()).size();++k)
260 {
261 // Nearest-neighbor interaction
262 static_cast<FieldPointBase<OtherParticleType,FieldType>* const>(lpt.bin(this->mpiRank())[k]) -> setZero();
263
264 typename SpatialCellType::CellMapType cellMap((lpt.bin(this->mpiRank())[k])->template neighborCells<_ParticleType>());
265
266 // loop over neighbor cells of the current particle
267 for (typename SpatialCellType::CellMapType::const_iterator cIter =cellMap.begin();
268 /* */ cIter!=cellMap.end();
269 /* */ ++cIter)
270 {
271 // loop over source particles in the neighbor cell
272 for(typename SpatialCellType::ParticleContainerType::const_iterator sIter=cIter->second->particleBegin();sIter!=cIter->second->particleEnd();++sIter) // loop over neighbor particles
273 {
274 // lpt.bin(this->mpiRank())[k]->template field<FieldType>() += FieldType::compute(**sIter,*lpt.bin(this->mpiRank())[k]);
275 if(static_cast<const SingleSourcePoint<ParticleType,FieldType>* const>(*sIter)->enabled)
276 {
277 *static_cast<FieldPointBase<OtherParticleType,FieldType>* const>(lpt.bin(this->mpiRank())[k]) += FieldType::compute(**sIter,*lpt.bin(this->mpiRank())[k]);
278 }
279 }
280 }
281
282 // Non nearest-neighbor interaction
283 if(FieldType::use_multipole)
284 {
285 typename SpatialCellType::CellMapType farCells(lpt.bin(this->mpiRank())[k]->template farCells<_ParticleType>());
286 // lpt.bin(this->mpiRank())[k]->template field<FieldType>() += FieldType::multipole(*lpt.bin(this->mpiRank())[k],farCells);
287 *static_cast<FieldPointBase<OtherParticleType,FieldType>* const>(lpt.bin(this->mpiRank())[k]) += FieldType::multipole(*lpt.bin(this->mpiRank())[k],farCells);
288
289 }
290
291 // Add contribution of other sources
292 if(sizeof...(OtherSourceTypes))
293 {
294 *static_cast<FieldPointBase<OtherParticleType,FieldType>* const>(lpt.bin(this->mpiRank())[k]) += FieldType::addSourceContribution(*lpt.bin(this->mpiRank())[k],otherSources...);
295 }
296
297 }
298
299
300 std::vector<int> interactionSizeVector(this->mpiProcs());
301 std::vector<int> interactionRankOffsetVector(this->mpiProcs());
302
303 size_t dofOffset(0);
304 for (int k=0;k<this->mpiProcs();++k)
305 {
306 // Store the number
307 interactionSizeVector[k]=lpt.bin(k).size()*FieldPointBase<OtherParticleType,FieldType>::DataPerParticle;
308 interactionRankOffsetVector[k]=dofOffset;
310 }
311
313 MPI_Allgatherv(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,
315 &interactionSizeVector[0],&interactionRankOffsetVector[0],
316 MPI_DOUBLE,MPI_COMM_WORLD);
317 }
318
319 /**********************************************************************/
320 template <typename OtherParticleType, typename FieldType,typename... OtherSourceTypes>
321 void computeField(OtherParticleType& part, const OtherSourceTypes&... otherSources) const
322 {
327 ParticleSystemBaseType::template computeField<OtherParticleType,FieldType,OtherSourceTypes...>(part,otherSources...);
328 }
329
330 /**********************************************************************/
331 template <class T>
332 friend T& operator<< (T& os, const ParticleContainerType& pS)
333 {
335 if (ModelMPIbase::mpiRank()==0)
336 {
337 for (typename ParticleContainerType::const_iterator pIter=pS.begin();pIter!=pS.end();++pIter)
338 {
339 os<<(*pIter);
340 }
341 }
342 return os;
343 }
344
345 };
346
347
348} // end namespace
349#endif
350
351
352// /**********************************************************************/
353// template <typename FieldType>
354// typename FieldType::ResultType getInteractionResult(const int& pID)
355// {
356// return FieldType::getResult(*(this->find(pID)->second));
357// }
358//
359// /**********************************************************************/
360// template <typename FieldType>
361// void resetInteraction()
362// {/*! Compute full particle interaction according to FieldType
363// */
364// for (typename ParticleContainerType::iterator iterJ=this->begin();iterJ!=this->end();++iterJ) // FOR SYMMETRIC INTERACTION iterJ STARTS AT iterI
365// {
366// FieldType::reset(*iterJ->second); // the constructor of FieldType actually computes the binary interaction between *iterI and *iterJ
367// }
368// }
369
370
371
372/**********************************************************************/
373// void partionByCells()
374// {
375// typedef typename SpatialCellObserverType::CellIdType CellIdType;
376// typedef std::map<CellIdType,const int,CompareVectorsByComponent<int,_ParticleType::dim> > ijk2idMapType;
377// // typedef std::map<Eigen::Matrix<int,_ParticleType::dim,1>,const int,CompareVectorsByComponent<int,_ParticleType::dim> > ijk2idMapType;
378// ijk2idMapType ijk2idMap;;
379// typedef std::map<int,const CellIdType> id2ijkMapType;
380// // typedef std::map<int,Eigen::Matrix<int,_ParticleType::dim,1> > id2ijkMapType;
381// id2ijkMapType id2ijkMap;;
382//
383// // Populate ijk2idMap and id2ijkMap
384// int k(0);
385// for (typename SpatialCellObserverType::CellMapType::const_iterator cIter =SpatialCellObserverType::cellBegin();
386// /* */ cIter!=SpatialCellObserverType::cellEnd();
387// /* */ ++cIter)
388// {
389// ijk2idMap.insert(std::pair<Eigen::Matrix<int,_ParticleType::dim,1>,const int>(cIter->second->cellID,k));
390// id2ijkMap.insert(std::pair<int,Eigen::Matrix<int,_ParticleType::dim,1> >(k,cIter->second->cellID));
391// ++k;
392// }
393//
394//
395//
396// // Compute weightScaleFactor to avoid int overflow in std::vector<int> vWeights
397// double totalWeightD(0.0);
398// for (typename SpatialCellObserverType::CellMapType::const_iterator cIter =SpatialCellObserverType::cellBegin();
399// /* */ cIter!=SpatialCellObserverType::cellEnd();
400// /* */ ++cIter)
401// {
402// totalWeightD+=cIter->second->n2WeightD();
403// }
404// double weightScaleFactor(1.0);
405// if (totalWeightD>=INT_MAX) // need to rescale actual "int" weights
406// {
407// weightScaleFactor=totalWeightD/INT_MAX;
408// if (this->mpiRank() == 0)
409// {
410//
411// std::cout<<"detected weight overflow. Rescaling weights: weightScaleFactor="<<weightScaleFactor<<std::endl;
412// }
413// }
414//
415// // Assemble std::vector<int> vWeights
416// std::vector<int> vWeights(SpatialCellObserverType::totalCells(),0); // vector containing the (scaled) weight of each cell
417// int kk(0);
418// int totalWeight(0);
419// for (typename SpatialCellObserverType::CellMapType::const_iterator cIter =SpatialCellObserverType::cellBegin();
420// /* */ cIter!=SpatialCellObserverType::cellEnd();
421// /* */ ++cIter)
422// {
423// vWeights[kk]=cIter->second->n2Weight()/weightScaleFactor;
424// totalWeight+=vWeights[kk];
425// ++kk;
426// }
427//
428// // Check overflow
429// if (totalWeightD>=INT_MAX) // actual totalWeight will overflow: need to rescale actual "int" weights
430// {
431// const double actualScaleFactor(totalWeightD/totalWeight);
432// if (this->mpiRank() == 0)
433// {
434// std::cout<<"totalWeight="<<totalWeight<<", totalWeightD="<<totalWeightD<<", actualScaleFactor="<<actualScaleFactor<<std::endl;
435// }
436// model_removeAssert(std::fabs(actualScaleFactor-weightScaleFactor)/weightScaleFactor<0.01 && "WEIGHT SCALE FACTOR MISMATCH");
437// }
438//
439//
440// // - Populate vector of neighbors and offsets
441// std::vector<int> vNeighbors;
442// std::vector<int> vOffsets;
443//
444// int offSet(0);
445// for (typename SpatialCellObserverType::CellMapType::const_iterator cIter =SpatialCellObserverType::cellBegin();
446// /* */ cIter!=SpatialCellObserverType::cellEnd();
447// /* */ ++cIter)
448// {
449// vOffsets.push_back(offSet);
450// for (typename SpatialCellObserverType::CellMapType::const_iterator nIter =cIter->second->neighborCellsBegin();
451// /* */ nIter!=cIter->second->neighborCellsEnd();
452// /* */ ++nIter)
453// {
454// if (cIter->second->cellID!=nIter->second->cellID)
455// {
456// typename ijk2idMapType::const_iterator mIter(ijk2idMap.find(nIter->second->cellID.transpose()));
457// vNeighbors.push_back(mIter->second);
458// ++offSet;
459// }
460// }
461// }
462// vOffsets.push_back(offSet); // push back an extra offset (required by Metis)
463//
464//
465// // Create a vector to store the output of the partitioner. Initialize at process=0
466// std::vector<int> partitionerRankVector(SpatialCellObserverType::totalCells(),0);
467//
468// // partition into the same number of MPI processes
469// int numPart = this->mpiProcs();
470//
471// if (numPart > 1)
472// {
473// int metisOptions[METIS_NOPTIONS];
474// int numVert = SpatialCellObserverType::totalCells(); // number of nodes
475// int ncon = 1; // number of weights per node
476// int objval;
477// float ubvec = 1.001; // load imbalance tolerance
478// METIS_SetDefaultOptions(metisOptions);
479// METIS_PartGraphKway(&numVert,&ncon,&vOffsets[0],&vNeighbors[0],
480// &vWeights[0],NULL,NULL,&numPart,NULL,
481// &ubvec,metisOptions,&objval,&partitionerRankVector[0]);
482// }
483//
484//
485//
486// // Assign each cell the corresponding processor, as found by metis
487// // int assignedParticleSize(0);
488// for (unsigned int i=0; i<SpatialCellObserverType::totalCells(); ++i)
489// {
490// typename id2ijkMapType::const_iterator idIter(id2ijkMap.find(i));
491// const CellIdType cellID(idIter->second);
492//
493// typename SpatialCellObserverType::isCellType isC(SpatialCellObserverType::isCell(cellID));
494// model_removeAssert(isC.first && "CELL MUST EXIST!");
495// isC.second->assignedRank=partitionerRankVector[i];
496//
497// if (this->mpiRank()==partitionerRankVector[i])
498// {
499// // assignedCells.push_back(isC.second);
500// // loop over particles in the current cell
501// for (typename SpatialCellType::ParticleContainerType::const_iterator pIter =isC.second->particleBegin();
502// /* */ pIter!=isC.second->particleEnd();
503// /* */ ++pIter)
504// {
505//
506// const typename ParticleContainerType::iterator nIter(this->find((*pIter)->sID));
507// // assert(nIter!=this->end());
508// assignedParticles.insert(std::make_pair((*pIter)->sID,nIter->second));
509// // assignedParticles.push_back(nIter->second);
510// }
511//
512//
513// // assignedParticleSize+=isC.second->size();
514// }
515// }
516//
517// // evaluate Performance
518// evaluatePerformance(partitionerRankVector,vWeights,totalWeight);
519//
520// // return assignedParticles.size();
521// // return assignedParticleSize;
522// }
523
524/**********************************************************************/
525// void partionByParticles()
526// {
527// // Initialize the vector containing the weight of each particle
528// std::vector<int> vWeights(this->size(),1);
529//
530// // Populate the vector containing the weight of each particle
531// int kk(0);
532// for (typename ParticleContainerType::const_iterator pIter =this->begin();
533// /* */ pIter!=this->end();
534// /* */ ++pIter)
535// {
536// // The weight of the kk-th particle is the n2Weight() of its cell
537// vWeights[kk]=pIter->second->pCell->n2Weight();
538// ++kk;
539// }
540//
541// // Initialize METIS vectors
542// std::vector<int> vNeighbors;
543// std::vector<int> vOffsets;
544//
545// // Populate METIS vectors
546// int offSet(0);
547// int pp(0);
548// for (typename ParticleContainerType::const_iterator pIter =this->begin();
549// /* */ pIter!=this->end();
550// /* */ ++pIter)
551// {
552// vOffsets.push_back(offSet);
553// vNeighbors.push_back(pp); // each particle is its only neighbor
554// ++offSet;
555// ++pp;
556// }
557// vOffsets.push_back(offSet); // push back an extra offset (required by Metis)
558//
559//
560// // partition into the same number of MPI processes
561// int numPart = this->mpiProcs();
562//
563// // Create a vector to store the output of the partitioner. Initialize at process=0
564// std::vector<int> partitionerRankVector(this->size(),0); //
565//
566// // Call METIS
567// if (numPart > 1)
568// {
569// int metisOptions[METIS_NOPTIONS];
570// int numVert = this->size(); // number of nodes
571// int ncon = 1; // number of weights per node
572// int objval;
573// float ubvec = 1.001; // load imbalance tolerance
574// METIS_SetDefaultOptions(metisOptions);
575// METIS_PartGraphKway(&numVert,&ncon,&vOffsets[0],&vNeighbors[0],
576// &vWeights[0],NULL,NULL,&numPart,NULL,
577// &ubvec,metisOptions,&objval,&partitionerRankVector[0]);
578// }
579//
580// // int assignedParticleSize(0);
581//
582// for (unsigned int i=0; i<this->size(); ++i)
583// {
584// if (this->mpiRank()==partitionerRankVector[i])
585// {
586// typename ParticleContainerType::iterator pIter(this->begin());
587// std::advance(pIter,i);
588// // assignedParticles.push_back(pIter->second);
589// assignedParticles.insert(std::make_pair(pIter->second->sID,pIter->second));
590//
591// // assignedParticleSize+=1;
592// }
593// }
594//
595// // return assignedParticleSize;
596// }
597
Template implementation of the Longest Processing Time (LPT) algorithm for task partiotioning.
void partition(const int &nProc)
static const int & mpiInitialized()
static const int & mpiRank()
ModelMPIbase(int argc, char *argv[])
static const int & mpiProcs()
SpatialCellObserver< _ParticleType, _ParticleType::dim > SpatialCellObserverType
SpatialCellObserverType::SpatialCellType SpatialCellType
_ParticleType::PositionType PositionType
ParticleType * addParticle(const PositionType &p, const AdditionalConstructorTypes &... args)
std::deque< _ParticleType > ParticleContainerType
ParticleSystemParallel(int argc, char *argv[])
void computeNeighborField(const OtherSourceTypes &... otherSources)
void computeField(OtherParticleType &part, const OtherSourceTypes &... otherSources) const
ParticleSystemBaseType::SpatialCellType SpatialCellType
ParticleSystemBaseType::ParticleContainerType ParticleContainerType
void computeField(std::deque< OtherParticleType > &fpDeq, const OtherSourceTypes &... otherSources) const
friend T & operator<<(T &os, const ParticleContainerType &pS)
ParticleType * addParticle(const typename ParticleSystemBaseType::PositionType &p, const AdditionalConstructorTypes &... args)
b
Definition int_R.m:16
A class template that contains the underlying FieldType data of the Derived particle.
static std::vector< Scalar > resultVector
static void resize(const unsigned int &k, const Scalar &val=Scalar())
int end[2]
void argv
void b int argc
vertex p[3]
int k