source: XIOS2/dev/mhedley/buildCompilationPatchesA/extern/remap/src/parallel_tree.cpp @ 2697

Last change on this file since 2697 was 2690, checked in by jderouillat, 7 weeks ago

Added the sorting of the localTree to ensure the reproductibilty of second order interpolation
File size: 17.1 KB
Line 
1#include <cassert>
2#include "node.hpp"
3#include "timerRemap.hpp"
4#include "circle.hpp"
5#include "meshutil.hpp"
6#include "polyg.hpp"
7#include "gridRemap.hpp"
8#include "intersect.hpp"
9#include "errhandle.hpp"
10#include "mpi_routing.hpp"
11#include "misc.hpp"
12
13#include "parallel_tree.hpp"
14
15namespace sphereRemap {
16
17static const int assignLevel = 2;
18
19// only the circle is packed, rest of node will be initialized on unpacking
20static void packNode(Node& node, char *buffer, int& index)
21{
22        if (buffer == NULL) // compute size only
23                index += 4 * sizeof(double);
24        else
25        {
26                *(Coord *)(&buffer[index]) = node.centre;
27                index += sizeof(Coord);
28
29                *(double *)(&buffer[index]) = node.radius;
30                index += sizeof(double);
31        }
32}
33
34static void unpackNode(Node& node, char *buffer, int& index)
35{
36        Coord centre;
37        double r;
38
39        centre = *(Coord *)(&buffer[index]);
40        index += sizeof(Coord);
41
42        r = *(double *)(&buffer[index]);
43        index += sizeof(double);
44
45        node.centre = centre;
46        node.radius = r;
47}
48
49
50static void packElement(Elt *ptElement, char *buffer, int& index)
51{
52        if (buffer == NULL) index += packedPolygonSize(*ptElement);
53        else packPolygon(*ptElement, buffer, index);
54}
55
56static void unpackElement(Elt *ptElement, char *buffer, int& index)
57{
58        unpackPolygon(*ptElement, buffer, index);
59}
60
61static void packVector(const vector<int>& vec, char *buf, int& pos)
62{
63        if (buf == NULL)
64                pos += sizeof(int) + vec.size()*sizeof(int);
65        else
66        {
67                *((int *) &(buf[pos])) = vec.size();
68                pos += sizeof(int);
69                for (int i = 0; i < vec.size(); i++)
70                {
71                        *((int *) &(buf[pos])) = vec[i];
72                        pos += sizeof(int);
73                }
74        }
75}
76
77static void unpackVector(vector<int>& vec, const char *buf, int& pos)
78{
79        vec.resize(*((int *) &(buf[pos])));
80        pos += sizeof(int);
81        for (int i = 0; i < vec.size(); i++)
82        {
83                vec[i] = *((int *) &(buf[pos]));
84                pos += sizeof(int);
85        }
86}
87
88static void assignRoute(CSampleTree& tree, const CCascadeLevel& cl)  // newroot <- root
89{
90        vector<int> routeRank(cl.group_size);
91        for (int i = 0; i < cl.group_size; i++)
92                routeRank[i] = i; //(cl.group_size + cl.polour() - i) % cl.group_size;
93        std::vector<int>::iterator rank = routeRank.begin();
94        tree.root->assignRoute(rank, assignLevel);
95}
96
97static void transferRoutedNodes(CMPIRouting& MPIRoute, /*const*/ vector<Node>& nodeSend, const vector<int>& route, vector<Node>& nodeRecv)
98{
99        /* `route` knows what goes where */
100        MPIRoute.init(route);
101        int nRecv = MPIRoute.getTotalSourceElement();
102        nodeRecv.resize(nRecv);
103        MPIRoute.transferToTarget(&nodeSend[0], &nodeRecv[0], packNode, unpackNode);
104}
105
106static void transferRoutedIntersections(CMPIRouting& MPIRoute, const vector<Node>& nodeSend, const vector<vector<int> >& route, vector<Node>& nodeRecv)
107{
108        // `route` knows what goes where
109        MPIRoute.init(route);
110        int nRecv = MPIRoute.getTotalSourceElement();
111        nodeRecv.resize(nRecv);
112        MPIRoute.transferToTarget((Node * /*mpi wants non-const*/) &nodeSend[0], &nodeRecv[0], packNode, unpackNode);
113//cout << MPIRoute.mpiRank << " ROUTE " << nRecv << ": " << nodeSend.size() << " " << nodeRecv.size() << "    ";
114}
115
116//CParallelTree::CParallelTree(MPI_Comm comm) : communicator(comm), cascade(MIN_NODE_SZ*MIN_NODE_SZ, comm)
117CParallelTree::CParallelTree(MPI_Comm comm) : communicator(comm), cascade(MAX_NODE_SZ*MAX_NODE_SZ*2, comm)
118{
119        treeCascade.reserve(cascade.num_levels);
120        for (int lev = 0; lev < cascade.num_levels; lev++)
121                treeCascade.push_back(CSampleTree(cascade.level[lev].group_size, assignLevel));
122}
123
124void CParallelTree::buildSampleTreeCascade(vector<Node>& sampleNodes /*route field will be modified*/, int level)
125{
126        buildSampleTree(treeCascade[level], sampleNodes, cascade.level[level]);
127        assignRoute(treeCascade[level] /*out*/, cascade.level[level] /*in*/);
128
129        if (level+1 < cascade.num_levels)
130        {
131                vector<int> route(sampleNodes.size());
132                treeCascade[level].routeNodes(route, sampleNodes, assignLevel);
133
134                vector<Node> routedNodes;
135                CMPIRouting mpiRoute(cascade.level[level].pg_comm);
136                transferRoutedNodes(mpiRoute, sampleNodes, route, routedNodes);
137                buildSampleTreeCascade(routedNodes, level+1);
138        }
139}
140
141void buildSampleTree(CSampleTree& tree, const vector<Node>& node, const CCascadeLevel& comm)
142/*
143        In the beginning all the sample elements are distributed
144        -> communicate to make available at each rank
145           so that each rank can build the same sample tree
146*/
147{
148        int n = node.size(); // number of samples intially on this rank (local)
149
150        int blocSize = comm.group_size * ipow(MAX_NODE_SZ, assignLevel);
151
152        int nrecv; // global number of samples  THIS WILL BE THE NUMBER OF LEAFS IN THE SAMPLE TREE
153        MPI_Allreduce(&n, &nrecv, 1, MPI_INT, MPI_SUM, comm.comm); // => size of sample tree does not depend on keepNodes!
154        double ratio = blocSize / (1.0 * nrecv);
155        int nsend = ratio * n + 1; // nsend = n_local_samples / n_global_samples * blocksize + 1 = blocksize/comm.size
156        if (nsend > n) nsend = n;
157
158        int *counts = new int[comm.size];
159        MPI_Allgather(&nsend, 1, MPI_INT, counts, 1, MPI_INT, comm.comm);
160
161        nrecv = 0;
162        int *displs = new int[comm.size];
163        for (int i = 0; i < comm.size; i++)
164        {
165                displs[i] = 4 * nrecv;
166                nrecv += counts[i];
167                counts[i] = 4 * counts[i];
168        }
169
170        /* pack circles around sample elements */
171        double *sendBuffer = new double[nsend*4];
172        int index = 0;
173        vector<int> randomArray(n);
174        randomizeArray(randomArray);
175        for (int i = 0; i < nsend; i++)
176        {
177                const Node& no = node[randomArray[i]];
178                *((Coord *) (sendBuffer + index)) = no.centre;
179                index += sizeof(Coord)/sizeof(*sendBuffer);
180                sendBuffer[index++] = no.radius;
181        }
182
183        /* each process needs the sample elements from all processes */
184        double *recvBuffer = new double[nrecv*4];
185        MPI_Allgatherv(sendBuffer, 4 * nsend, MPI_DOUBLE, recvBuffer, counts, displs, MPI_DOUBLE, comm.comm);
186        delete[] sendBuffer;
187        delete[] counts;
188        delete[] displs;
189
190        /* unpack */
191/*
192        randomArray.resize(nrecv);
193        randomizeArray(randomArray);
194        tree.leafs.resize(nrecv);
195        index = 0;
196
197
198  for (int i = 0; i < nrecv; i++)
199        {
200                Coord x = *(Coord *)(&recvBuffer[index]);
201                index += sizeof(Coord)/sizeof(*recvBuffer);
202                double radius = recvBuffer[index++];
203                tree.leafs[randomArray[i]].centre = x;
204                tree.leafs[randomArray[i]].radius = radius;
205
206        }
207*/
208
209  randomArray.resize(blocSize);
210        randomizeArray(randomArray);
211        tree.leafs.resize(blocSize);
212        index = 0;
213 
214  size_t s=(sizeof(Coord)/sizeof(*recvBuffer)+1)*nrecv ;
215 
216  for (int i = 0; i < blocSize; i++)
217        {
218                Coord x = *(Coord *)(&recvBuffer[index%s]);
219                index += sizeof(Coord)/sizeof(*recvBuffer);
220                double radius = recvBuffer[index%s];
221    index++ ;
222                tree.leafs[randomArray[i]].centre = x;
223                tree.leafs[randomArray[i]].radius = radius;
224
225        }
226
227
228        delete [] recvBuffer;
229
230        CTimer::get("buildSampleTree(local)").resume();
231        tree.build(tree.leafs);
232//      cout << "SampleTree build : assign Level " << assignLevel << " nb Nodes : " << tree.levelSize[assignLevel] << endl;
233        CTimer::get("buildSampleTree(local)").suspend();
234        CTimer::get("buildSampleTree(local)").print();
235
236        /* End gracefully if sample tree could not be constructed with desired number of nodes on assignLevel */
237        int allok, ok = (tree.levelSize[assignLevel] == comm.group_size);
238        if (!ok)
239  {
240    cerr << comm.rank << ": PROBLEM: (node assign)" << tree.levelSize[assignLevel] << " != " << comm.group_size << " (keepNodes)" 
241         << "   node size : "<<node.size()<<"   bloc size : "<<blocSize<<"  total number of leaf : "<<tree.leafs.size()<<endl ;
242/*
243        MPI_Allreduce(&ok, &allok, 1, MPI_INT, MPI_PROD, communicator);
244        if (!allok) {
245                MPI_Finalize();
246                exit(1);
247        }
248*/
249    MPI_Abort(MPI_COMM_WORLD,-1) ;
250  }
251/*
252  cout<<"*******************************************"<<endl ;
253  cout<<"******* Sample Tree output        *********"<<endl ;
254  cout<<"*******************************************"<<endl ;
255  tree.output(cout,1) ;
256
257  cout<<"*******************************************"<<endl ;
258*/
259
260  assert(tree.root->incluCheck() == 0);
261}
262
263struct compareElts
264{
265  bool operator()(const std::tuple<double, double, double>& lhs, const std::tuple<double, double, double>& rhs) const
266  {
267    if (get<0>(lhs) < get<0>(rhs)) return true;
268    else if (get<0>(lhs) > get<0>(rhs)) return false;
269    else {
270      if (get<1>(lhs) < get<1>(rhs)) return true;
271      else if (get<1>(lhs) > get<1>(rhs)) return false;
272      else {
273        if (get<2>(lhs) < get<2>(rhs)) return true;
274        else if (get<2>(lhs) > get<2>(rhs)) return false;
275        else return true;
276      }
277    }
278  }
279};
280
281void CParallelTree::buildLocalTree(const vector<Node>& node, const vector<int>& route)
282{
283        CMPIRouting MPIRoute(communicator);
284        MPI_Barrier(communicator);
285        CTimer::get("buildLocalTree(initRoute)").resume();
286        MPIRoute.init(route);
287        CTimer::get("buildLocalTree(initRoute)").suspend();
288        CTimer::get("buildLocalTree(initRoute)").print();
289
290        nbLocalElements = MPIRoute.getTotalSourceElement();
291        localElements = new Elt[nbLocalElements];
292
293        vector<Elt*> ptElement(node.size());
294        for (int i = 0; i < node.size(); i++)
295                ptElement[i] = (Elt *) (node[i].data);
296
297        vector<Elt*> ptLocalElement(nbLocalElements);
298        for (int i = 0; i < nbLocalElements; i++)
299                ptLocalElement[i] = &localElements[i];
300
301       
302        CTimer::get("buildLocalTree(transfer)").resume();
303        MPIRoute.transferToTarget(&ptElement[0], &ptLocalElement[0], packElement, unpackElement);
304        CTimer::get("buildLocalTree(transfer)").suspend();
305        CTimer::get("buildLocalTree(transfer)").print();
306
307        CTimer::get("buildLocalTree(local)").resume();
308       
309        // Force order access to localElements for reproductibilty using the map below
310        std::map<std::tuple<double, double, double>, int, compareElts> indexesFromCoords;
311        for (int i = 0; i < nbLocalElements; i++)
312        {
313            Elt& elt = localElements[i];
314            indexesFromCoords[{elt.x.x,elt.x.y,elt.x.z}] = i;
315        }
316
317        int mpiRank;
318        MPI_Comm_rank(communicator, &mpiRank);
319        localTree.leafs.reserve(nbLocalElements);
320        for (const auto& itIndex : indexesFromCoords)
321        {
322                int i = itIndex.second; 
323                Elt& elt = localElements[i];
324                elt.id.ind = i;
325                elt.id.rank = mpiRank;
326                localTree.leafs.push_back(Node(elt.x, cptRadius(elt), &localElements[i]));
327        }
328        localTree.build(localTree.leafs);
329
330        cptAllEltsGeom(localElements, nbLocalElements, srcGrid.pole);
331        CTimer::get("buildLocalTree(local)").suspend();
332        CTimer::get("buildLocalTree(local)").print();
333}
334
335void CParallelTree::build(vector<Node>& node, vector<Node>& node2)
336{
337
338        int assignLevel = 2;
339        int nbSampleNodes = 2*ipow(MAX_NODE_SZ + 1, assignLevel);
340
341
342  long int nb1, nb2, nb, nbTot ;
343  nb1=node.size() ; nb2=node2.size() ;
344  nb=nb1+nb2 ;
345  MPI_Allreduce(&nb, &nbTot, 1, MPI_LONG, MPI_SUM, communicator) ;
346  int commSize ;
347  MPI_Comm_size(communicator,&commSize) ;
348 
349        // make multiple of two
350        nbSampleNodes /= 2;
351        nbSampleNodes *= 2;
352//  assert( nbTot > nbSampleNodes*commSize) ;
353   
354  int nbSampleNodes1 = nbSampleNodes * (nb1*commSize)/(1.*nbTot) ;
355  int nbSampleNodes2 = nbSampleNodes * (nb2*commSize)/(1.*nbTot) ;
356 
357
358//      assert(node.size() > nbSampleNodes);
359//      assert(node2.size() > nbSampleNodes);
360//      assert(node.size() + node2.size() > nbSampleNodes);
361        vector<Node> sampleNodes; sampleNodes.reserve(nbSampleNodes1+nbSampleNodes2);
362
363        vector<int> randomArray1(node.size());
364        randomizeArray(randomArray1);
365        vector<int> randomArray2(node2.size());
366        randomizeArray(randomArray2);
367
368/*     
369        int s1,s2 ;
370       
371        if (node.size()< nbSampleNodes/2)
372        {
373          s1 = node.size() ;
374          s2 = nbSampleNodes-s1 ;
375        }
376        else if (node2.size()< nbSampleNodes/2)
377        {
378          s2 = node.size() ;
379          s1 = nbSampleNodes-s2 ;
380        }
381        else
382        {
383          s1=nbSampleNodes/2 ;
384          s2=nbSampleNodes/2 ;
385        }
386*/
387        for (int i = 0; i <nbSampleNodes1; i++) sampleNodes.push_back(Node(node[randomArray1[i%nb1]].centre,  node[randomArray1[i%nb1]].radius, NULL));
388        for (int i = 0; i <nbSampleNodes2; i++) sampleNodes.push_back(Node(node2[randomArray2[i%nb2]].centre, node2[randomArray2[i%nb2]].radius, NULL));
389
390/*         
391        for (int i = 0; i < nbSampleNodes/2; i++)
392        {
393          sampleNodes.push_back(Node(node[randomArray1[i]].centre,  node[randomArray1[i]].radius, NULL));
394          sampleNodes.push_back(Node(node2[randomArray2[i]].centre, node2[randomArray2[i]].radius, NULL));
395        }
396*/
397        CTimer::get("buildParallelSampleTree").resume();
398        //sampleTree.buildParallelSampleTree(sampleNodes, cascade);
399        buildSampleTreeCascade(sampleNodes);
400        CTimer::get("buildParallelSampleTree").suspend();
401        CTimer::get("buildParallelSampleTree").print();
402
403        //route source mesh
404        CTimer::get("parallelRouteNode").resume();
405        vector<int> route(node.size());
406        routeNodes(route /*out*/, node);
407        CTimer::get("parallelRouteNode").suspend();
408        CTimer::get("parallelRouteNode").print();
409
410        CTimer::get("buildLocalTree").resume();
411        buildLocalTree(node, route);
412        CTimer::get("buildLocalTree").suspend();
413        CTimer::get("buildLocalTree").print();
414
415        CTimer::get("buildRouteTree").resume();
416        /* update circles of tree cascade so it can be used for routing */
417        updateCirclesForRouting(localTree.root->centre, localTree.root->radius);
418        CTimer::get("buildRouteTree").suspend();
419        CTimer::get("buildRouteTree").print();
420}
421
422void CParallelTree::routeNodes(vector<int>& route, vector<Node>& nodes /*route field used*/, int level)
423{
424        treeCascade[level].routeNodes(route /*assign*/, nodes, assignLevel);
425
426        if (level+1 < cascade.num_levels)
427        {
428                vector<Node> routedNodes;
429                CMPIRouting MPIRoute(cascade.level[level].pg_comm);
430                transferRoutedNodes(MPIRoute, nodes, route /*use*/, routedNodes);
431                vector<int> globalRank(routedNodes.size());
432                routeNodes(globalRank, routedNodes, level + 1);
433                MPIRoute.transferFromSource(&route[0] /*override*/, &globalRank[0]);
434        }
435        else
436        {
437                CMPIRouting MPIRoute(cascade.level[level].comm); // or use pg_comm, on last cascade level identical
438                MPIRoute.init(route);
439                int nbRecvNode = MPIRoute.getTotalSourceElement();
440                vector<int> globalRank(nbRecvNode);
441                for (int i = 0; i < globalRank.size(); i++)
442                        globalRank[i] = cascade.level[0].rank;
443                MPIRoute.transferFromSource(&route[0] /*override*/, &globalRank[0]);
444        }
445}
446
447/* assume `to` to be empty vector at entry */
448void linearize(const vector<vector<int> >& from, vector<int>& to)
449{
450        int cnt = 0;
451        for (int i = 0; i < from.size(); ++i)
452                cnt += from[i].size();
453        to.resize(cnt);
454        vector<int>::iterator dest = to.begin();
455        for (int i = 0; i < from.size(); ++i)
456                dest = copy(from[i].begin(), from[i].end(), dest);
457}
458
459/* at entry `to` must already have it's final shape and only values are overriden */
460void delinearize(const vector<int>& from, vector<vector<int> >& to)
461{
462        vector<int>::const_iterator end, src = from.begin();
463        for (int i = 0; i < to.size(); ++i, src=end)
464                copy(src, end = src + to[i].size(), to[i].begin());
465}
466
467void CParallelTree::routeIntersections(vector<vector<int> >& routes, vector<Node>& nodes, int level)
468{
469        treeCascade[level].routeIntersections(routes /*assign*/, nodes);
470
471        if (level+1 < cascade.num_levels)
472        {
473                vector<Node> routedNodes;
474                CMPIRouting MPIRoute(cascade.level[level].pg_comm);
475
476                vector<int> flattenedRoutes1;
477                linearize(routes, flattenedRoutes1);
478                vector<Node> double_nodes(flattenedRoutes1.size());
479                int j = 0;
480                for (int i = 0; i < routes.size(); ++i)
481                        for (int k = 0; k < routes[i].size(); ++k, ++j)
482                                double_nodes[j] = nodes[i];
483                transferRoutedNodes(MPIRoute, double_nodes, flattenedRoutes1 /*use*/, routedNodes);
484                vector<vector<int> > globalRanks(routedNodes.size());
485                routeIntersections(globalRanks /*out*/, routedNodes /*in*/, level + 1);
486                vector<vector<int> > flattenedRoutes(flattenedRoutes1.size());
487                // transferFromSource expects sizes (nbSendNode=flattenedRoutes, nbRecvNode=routedNodes.size())
488                MPIRoute.transferFromSource(&flattenedRoutes[0], &globalRanks[0], packVector, unpackVector);
489                for (int i = 0, j = 0; i < routes.size(); ++i)
490                {
491                        int old_size = routes[i].size();
492                        routes[i].resize(0);
493                        for (int k = 0; k < old_size; ++k, ++j)
494                                for (int l = 0; l < flattenedRoutes[j].size(); ++l)
495                                        routes[i].push_back(flattenedRoutes[j][l]);
496                }
497                assert(j == flattenedRoutes1.size());
498
499        }
500        else
501        {
502                CMPIRouting MPIRoute(cascade.level[level].comm);
503                MPIRoute.init(routes);
504                int nbRecvNode = MPIRoute.getTotalSourceElement();
505                vector<int> globalRanks(nbRecvNode, cascade.level[0].rank);
506                vector<int> flattened_routes;
507                linearize(routes, flattened_routes);
508                MPIRoute.transferFromSource(&flattened_routes[0], &globalRanks[0]);
509                delinearize(flattened_routes, routes);
510        }
511        if (level!=level)
512        {
513                for (int i = 0; i < routes.size(); ++i)
514                        for (int k = 0; k < routes[i].size(); ++k)
515                                if (routes[i][k] == cascade.level[0].rank) routes[i].erase(routes[i].begin() + (k--));
516        }
517}
518
519void CParallelTree::updateCirclesForRouting(Coord rootCentre, double rootRadius, int level)
520{
521        if (level + 1 < cascade.num_levels) // children in the cascade have to update first
522        {
523                updateCirclesForRouting(rootCentre, rootRadius, level + 1);
524                rootCentre = treeCascade[level+1].root->centre;
525                rootRadius = treeCascade[level+1].root->radius;
526        }
527
528        // gather circles on this level of the cascade
529        int pg_size;
530        MPI_Comm_size(cascade.level[level].pg_comm, &pg_size);
531        vector<Coord> allRootCentres(pg_size);
532        vector<double> allRootRadia(pg_size);
533        MPI_Allgather(&rootCentre, 3, MPI_DOUBLE, &allRootCentres[0], 3, MPI_DOUBLE, cascade.level[level].pg_comm);
534        MPI_Allgather(&rootRadius, 1, MPI_DOUBLE, &allRootRadia[0],   1, MPI_DOUBLE, cascade.level[level].pg_comm);
535
536        // now allRootsRadia and allRootCentres must be inserted into second levels of us and propagated to root
537        treeCascade[level].root->assignCircleAndPropagateUp(&allRootCentres[0], &allRootRadia[0], assignLevel);
538}
539
540CParallelTree::~CParallelTree()
541{
542        delete [] localElements;
543}
544
545}
Note: See TracBrowser for help on using the repository browser.