← Back to branch summary

~fluidity-core/fluidity/embedded_models

« back to all changes in this revision

Viewing changes to libadaptivity/load_balance/src/formHalo2.cpp

  • Committer: Timothy Bond
  • Date: 2011-04-14 15:40:24 UTC
  • Revision ID: timothy.bond@imperial.ac.uk-20110414154024-116ci9gq6mwigmaw
Following the move from svn to bzr we change the nature of inclusion of these
four software libraries. Previously, they were included as svn externals and
pulled at latest version for trunk, pinned to specific versions for release
and stable trunk. Since bzr is less elegant at dealing with externals we have
made the decision to include the packages directly into the trunk instead.

At this import the versions are:

libadaptivity: r163
libvtkfortran: r67
libspud: r545
libmba2d: r28

Show diffs side-by-side

added added

removed removed

Lines of Context:
libadaptivity/load_balance/src/formHalo2.cpp
 
1
/* Copyright (C) 2006 Imperial College London and others.
 
2
 
 
3
 Please see the AUTHORS file in the main source directory for a full list
 
4
 of copyright holders.
 
5
 
 
6
 Dr Gerard J Gorman
 
7
 Applied Modelling and Computation Group
 
8
 Department of Earth Science and Engineering
 
9
 Imperial College London
 
10
 
 
11
 g.gorman@imperial.ac.uk
 
12
 
 
13
 This library is free software; you can redistribute it and/or
 
14
 modify it under the terms of the GNU Lesser General Public
 
15
 License as published by the Free Software Foundation; either
 
16
 version 2.1 of the License.
 
17
 
 
18
 This library is distributed in the hope that it will be useful,
 
19
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 
20
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
21
 Lesser General Public License for more details.
 
22
 
 
23
 You should have received a copy of the GNU Lesser General Public
 
24
 License along with this library; if not, write to the Free Software
 
25
 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
 
26
 USA
 
27
*/
 
28
#include <mpi.h>
 
29
#include <assert.h>
 
30
#include <vector>
 
31
#include <deque>
 
32
#include <string>
 
33
#include <map>
 
34
#include <set>
 
35
 
 
36
using namespace std;
 
37
 
 
38
#include "Mesh.h"
 
39
#include "packing.h"
 
40
#include "c++debug.h"
 
41
#include "comTools.h"
 
42
 
 
43
//#define CREATE_VELOCITY_HALO2 1
 
44
 
 
45
void Mesh::formHalo2(){
 
46
#ifndef NDEBUG
 
47
  unsigned len45 = element_list.size();
 
48
  for(unsigned i=0; i<len45; i++)
 
49
    for(unsigned j=0; j<len45; j++)
 
50
      if(i!=j)
 
51
        assert(element_list[i] != element_list[j]);
 
52
  #endif
 
53
 
 
54
  // This is what we are looking for.
 
55
  vector< deque<unsigned> > halo2Elements(NProcs);
 
56
  set<unsigned> dangerElements;
 
57
  {  
 
58
    // Wizz through all elements.
 
59
    unsigned pos = 0;
 
60
    for(deque<Element>::const_iterator elm = element_list.begin(); elm != element_list.end(); ++elm, ++pos){
 
61
 
 
62
      const vector<unn_t>& nodes = (*elm).get_enlist();
 
63
      vector<bool> halo1(NProcs, false);
 
64
      bool interesting=false;
 
65
 
 
66
      for(int p=0; p<NProcs; p++){
 
67
        if(p==MyRank) continue;
 
68
        
 
69
        // Record if elm has a halo node with p.
 
70
        for( vector<unn_t>::const_iterator jt = nodes.begin(); jt!=nodes.end(); ++jt){
 
71
          halo1[p] = (halo_nodes[p].find(*jt)!= halo_nodes[p].end());
 
72
          
 
73
          // Early exit...
 
74
          if(halo1[p])
 
75
            break;
 
76
        }
 
77
        
 
78
        if(!halo1[p]){
 
79
          // Record if elm shares a node with p, but is not itself a halo1 element for p.
 
80
          for(vector<unn_t>::const_iterator jt=nodes.begin(); jt!=nodes.end(); ++jt){
 
81
            if(shared_nodes[p].find(*jt)!=shared_nodes[p].end()){
 
82
              interesting = true;
 
83
              halo2Elements[p].push_back(pos);
 
84
              break;
 
85
            }
 
86
          }
 
87
        }
 
88
 
 
89
      }
 
90
      
 
91
      if(interesting){
 
92
        // A halo2 element may be multiably sent if it contains halo1
 
93
        // nodes for any domain.
 
94
        for(unsigned p=0;p<(unsigned)NProcs;p++){
 
95
          if(halo1[p]){
 
96
            dangerElements.insert( pos );
 
97
            break;
 
98
          }
 
99
        }         
 
100
      }
 
101
    }
 
102
  } // Have halo2 elements and a set of dangerElements.
 
103
  
 
104
 
 
105
  //
 
106
  // Identify all the nodes that must be sent.
 
107
  //
 
108
  vector< set<unsigned> > sendhalo2nodes(NProcs);
 
109
  vector< set<unsigned> > sendhalo2pnodes(NProcs);
 
110
  for(unsigned p=0; p<(unsigned)NProcs; p++){
 
111
 
 
112
    for(deque<unsigned>::const_iterator elm=halo2Elements[p].begin(); elm!=halo2Elements[p].end(); ++elm){
 
113
      // It's enough to just check the volume elements
 
114
      if(element_list[*elm].get_flags() & ELM_SURFACE)
 
115
        continue;
 
116
      
 
117
      {// Regular nodes 
 
118
        const vector<unn_t>& nodes = element_list[*elm].get_enlist();
 
119
        for(vector<unn_t>::const_iterator nod=nodes.begin(); nod!=nodes.end(); ++nod){
 
120
          sendhalo2nodes[p].insert( *nod );
 
121
        }
 
122
      }
 
123
      {// Pressure nodes
 
124
        const vector<unn_t>& nodes = element_list[*elm].get_MFenlist();
 
125
        for(vector<unn_t>::const_iterator nod=nodes.begin(); nod!=nodes.end(); ++nod){
 
126
          sendhalo2pnodes[p].insert( *nod );
 
127
        }
 
128
      }
 
129
    }
 
130
    
 
131
    {// Remove nodes that p should already know through halo1.
 
132
      set<unsigned> toDel;
 
133
      for(set<unsigned>::const_iterator it=sendhalo2nodes[p].begin(); it != sendhalo2nodes[p].end(); ++it){
 
134
        if(shared_nodes[p].find( *it ) != shared_nodes[p].end()){
 
135
          toDel.insert( *it );
 
136
        }
 
137
      }
 
138
      for(set<unsigned>::const_iterator it=toDel.begin(); it != toDel.end(); ++it){
 
139
        sendhalo2nodes[p].erase( *it );
 
140
      }
 
141
    }
 
142
    
 
143
    {// Remove pressure nodes that p should already know through halo1.
 
144
      set<unsigned> toDel;
 
145
      for(set<unsigned>::const_iterator it=sendhalo2pnodes[p].begin(); it != sendhalo2pnodes[p].end(); ++it){
 
146
        if( shared_pnodes[p].find( *it ) != shared_pnodes[p].end() ){
 
147
          toDel.insert( *it );
 
148
        }
 
149
      }   
 
150
      for(set<unsigned>::const_iterator it=toDel.begin(); it != toDel.end(); ++it){
 
151
        sendhalo2pnodes[p].erase( *it );
 
152
      }
 
153
    }
 
154
  }
 
155
 
 
156
  //
 
157
  // At this point we have identified all the information which we
 
158
  // want to communicate: 
 
159
  // vector< deque<unsigned> > elems2send( NProcs );
 
160
  // vector< set<unsigned> > nodes2send( NProcs );
 
161
  //
 
162
  
 
163
  // Make the send-packs
 
164
  vector< vector<char> > SendRecvBuffer(NProcs);
 
165
  
 
166
  { // Allocate space for buffers
 
167
    unsigned max_node_size      = max_nodepack_size();
 
168
    unsigned max_pnode_size     = max_pressurepack_size();
 
169
    unsigned max_element_size   = max_elementpack_size();
 
170
    unsigned space_for_unsigned = MPI::UNSIGNED.Pack_size(1, MPI::COMM_WORLD);
 
171
    
 
172
    for(int i=0; i<NProcs; i++){
 
173
      unsigned nbytes = space_for_unsigned          +
 
174
        space_for_unsigned*halo2Elements[i].size() +
 
175
        space_for_unsigned                          +
 
176
        max_element_size*halo2Elements[i].size()    +
 
177
        space_for_unsigned                          +
 
178
        max_node_size*sendhalo2nodes[i].size()      +
 
179
        space_for_unsigned                          +
 
180
        max_pnode_size*sendhalo2pnodes[i].size();
 
181
      
 
182
      SendRecvBuffer[i].resize( (unsigned)(1.1*nbytes) );
 
183
    }
 
184
  }
 
185
  
 
186
  vector<int> offsets(NProcs, 0); // int because of mpi calls
 
187
  
 
188
  // Pack.
 
189
  for(int i=0; i<NProcs; i++){
 
190
    int len = SendRecvBuffer[i].size();
 
191
    
 
192
    if( (i == MyRank)||(len == 0) )
 
193
      continue;
 
194
    
 
195
    char *buff = &(SendRecvBuffer[i][0]);
 
196
    
 
197
    // Elements
 
198
    unsigned cnt = halo2Elements[i].size();
 
199
    MPI::UNSIGNED.Pack(&cnt, 1, buff, len, offsets[i], MPI::COMM_WORLD);
 
200
    
 
201
    ECHO("Packing "<<cnt<<" halo2 elements for "<<i<<".");
 
202
    for(unsigned j=0; j<cnt; j++){
 
203
 
 
204
      // Dangerious?
 
205
      unsigned danger = 0;
 
206
      if( dangerElements.find( halo2Elements[i][j] ) != dangerElements.end() )
 
207
        danger = 1;
 
208
      MPI::UNSIGNED.Pack(&danger, 1, buff, len, offsets[i], MPI::COMM_WORLD);
 
209
 
 
210
      // Pack element
 
211
      element_list[ halo2Elements[i][j] ].pack(buff, len, offsets[i]);
 
212
    }
 
213
    // Nodes
 
214
    cnt = sendhalo2nodes[i].size();
 
215
    MPI::UNSIGNED.Pack(&cnt, 1, buff, len, offsets[i], MPI::COMM_WORLD);
 
216
 
 
217
    ECHO("Packing "<<cnt<<" halo2 nodes for "<<i<<".");    
 
218
    for(set<unsigned>::const_iterator it=sendhalo2nodes[i].begin(); it!=sendhalo2nodes[i].end(); ++it){
 
219
      const Node& node = node_list.unn( *it );
 
220
      node.pack(buff, len, offsets[i]);
 
221
    }
 
222
    
 
223
    // Pressure nodes
 
224
    cnt = sendhalo2pnodes[i].size();
 
225
    MPI::UNSIGNED.Pack(&cnt, 1, buff, len, offsets[i], MPI::COMM_WORLD);
 
226
    ECHO("Packing "<<cnt<<" halo2 pressure nodes for "<<i<<".");
 
227
 
 
228
    for(set<unsigned>::const_iterator it=sendhalo2pnodes[i].begin(); it!=sendhalo2pnodes[i].end(); ++it){
 
229
      const PressureNode& node = MFnode_list.unn( *it );
 
230
      node.pack(buff, len, offsets[i]);
 
231
    }
 
232
    
 
233
    assert(offsets[i] <= (int)SendRecvBuffer[i].size());
 
234
  }
 
235
  
 
236
  // Clean-up.
 
237
  dangerElements.clear();
 
238
  halo2Elements.clear();
 
239
  sendhalo2nodes.clear();
 
240
  sendhalo2pnodes.clear();
 
241
  
 
242
  // Send/recieve everything.
 
243
  allSendRecv(SendRecvBuffer);
 
244
 
 
245
  { // Unpacking.
 
246
    ECHO("Starting unpacking.");
 
247
    
 
248
    for(int p=0; p<NProcs; p++){
 
249
      offsets[p] = 0;
 
250
    }
 
251
 
 
252
    set<Element>  DangerElements;
 
253
    set<unsigned> ReceivedNodes;
 
254
    set<unsigned> ReceivedPressureNodes;
 
255
    vector< set<unsigned> > extraHalo(NProcs);
 
256
    
 
257
    // Paranoid
 
258
    assert(SendRecvBuffer.size() == (unsigned)NProcs);
 
259
 
 
260
#ifndef NDEBUG
 
261
    CHECK( element_list.size() );       
 
262
    do_element_headcount();
 
263
    CHECK( num_elements("total") );
 
264
    CHECK( num_elements("volume") );
 
265
    CHECK( num_elements("surface") );
 
266
#endif
 
267
   
 
268
    for(int p=0; p<NProcs; p++){
 
269
      int nbytes = SendRecvBuffer[p].size();
 
270
      if( (p == MyRank)||(nbytes == 0) )
 
271
        continue;
 
272
      
 
273
      char *buffer = &(SendRecvBuffer[p][0]);
 
274
      { // elements
 
275
        unsigned cnt;
 
276
        MPI::UNSIGNED.Unpack(buffer, nbytes, &cnt, 1, offsets[p], MPI::COMM_WORLD);
 
277
        
 
278
        ECHO("Unpacking "<<cnt<<" elements from "<<p<<".");
 
279
        
 
280
        for(unsigned j=0; j<cnt; j++){
 
281
          Element element;
 
282
          
 
283
          ECHO("Unpacking "<<j<<"...");
 
284
          
 
285
          // Unpack danger flag.
 
286
          unsigned danger;
 
287
          MPI::UNSIGNED.Unpack(buffer, nbytes, &danger, 1, offsets[p], MPI::COMM_WORLD);
 
288
          CHECK(danger);
 
289
          
 
290
          // Unpack element...      
 
291
          element.unpack(buffer, nbytes, offsets[p]);
 
292
          ECHO("unpacked.");
 
293
          
 
294
          if(danger){
 
295
            ECHO("Danger element...taking evasive action."); 
 
296
            DangerElements.insert( element );         
 
297
          }else{
 
298
            element_list.push_back(element);
 
299
          }  
 
300
        }
 
301
      } // finished unpacking elements.
 
302
      
 
303
      { // nodes
 
304
        unsigned cnt;
 
305
        MPI::UNSIGNED.Unpack(buffer, nbytes, &cnt, 1, offsets[p], MPI::COMM_WORLD);
 
306
        ECHO("Unpacking "<<cnt<<" nodes from "<<p<<".");
 
307
        
 
308
        for(unsigned j=0; j<cnt; j++){
 
309
          Node node;
 
310
          node.unpack(buffer, nbytes, offsets[p]);
 
311
          unsigned unn   = node.get_unn();
 
312
          unsigned owner = node.get_owner();
 
313
 
 
314
          if(halo_nodes[owner].find( unn ) == halo_nodes[owner].end())
 
315
            if(ReceivedNodes.find( unn ) == ReceivedNodes.end()){
 
316
              assert( !node_list.contains_unn(unn) );
 
317
              node_list.push_back( node );
 
318
              ReceivedNodes.insert( unn );
 
319
            }
 
320
        }
 
321
        
 
322
      } // finished unpacking nodes
 
323
      
 
324
      { // pressure nodes
 
325
        unsigned cnt;
 
326
        MPI::UNSIGNED.Unpack(buffer, nbytes, &cnt, 1, offsets[p], MPI::COMM_WORLD);
 
327
        ECHO("Unpacking "<<cnt<<" pressure nodes from "<<p<<".");
 
328
        
 
329
        for(unsigned j=0; j<cnt; j++){
 
330
          PressureNode node;
 
331
          node.unpack(buffer, nbytes, offsets[p]);
 
332
          unsigned unn = node.get_unn();
 
333
          unsigned owner = node.get_owner();
 
334
 
 
335
          if(halo_pnodes[owner].find( unn ) == halo_pnodes[owner].end())
 
336
            if(ReceivedPressureNodes.find( node.get_unn() ) == ReceivedPressureNodes.end()){
 
337
              assert( !MFnode_list.contains_unn(unn) );
 
338
              MFnode_list.push_back( node );
 
339
              ReceivedPressureNodes.insert( node.get_unn() );
 
340
              
 
341
              // This next line is tricky as "node.get_owner()" is used
 
342
              // rather than p. This is because the owner of this node is
 
343
              // not necessarly the same as the processor that sent
 
344
              // it.
 
345
              extraHalo[owner].insert( node.get_unn() );
 
346
            }
 
347
        }
 
348
        
 
349
      } // finished unpacking pressure nodes
 
350
      
 
351
    }
 
352
  
 
353
    //
 
354
    // Finally, add in the danger elements
 
355
    //
 
356
    if( !DangerElements.empty() ){
 
357
      ECHO("Add in danger elements.");
 
358
      
 
359
      for(set<Element>::const_iterator it=DangerElements.begin(); it!=DangerElements.end(); ++it)
 
360
        element_list.push_back( *it ); 
 
361
      DangerElements.clear();
 
362
    }
 
363
    
 
364
#ifndef NDEBUG
 
365
    CHECK( element_list.size() );       
 
366
    do_element_headcount();
 
367
    CHECK( num_elements("total") );
 
368
    CHECK( num_elements("volume") );
 
369
    CHECK( num_elements("surface") );
 
370
#endif
 
371
 
 
372
    //
 
373
    // Communicate extra halo stuff.
 
374
    //
 
375
    vector< set<unsigned> > extraShared(NProcs);
 
376
    ECHO("Communicating the extended halo...");
 
377
    assert(extraHalo[MyRank].size() == 0);    
 
378
    allSendRecv(extraHalo, extraShared);
 
379
    assert(extraShared[MyRank].size() == 0);
 
380
    ECHO("...done.");
 
381
    
 
382
    // 
 
383
    // Add in extra halo values...
 
384
    //
 
385
    for(int p=0; p<NProcs; p++){
 
386
      if(p==MyRank)
 
387
        continue;
 
388
      
 
389
      // add in extra halo values
 
390
      ECHO("Adding halo nodes for "<<p);
 
391
      halo_pnodes[p].insert(extraHalo[p].begin(), extraHalo[p].end());
 
392
      extraHalo[p].clear();
 
393
      
 
394
      ECHO("Adding shared nodes for "<<p);
 
395
      shared_pnodes[p].insert(extraShared[p].begin(), extraShared[p].end());
 
396
      extraShared[p].clear();
 
397
      
 
398
    }
 
399
    extraHalo.clear();
 
400
    extraShared.clear();
 
401
    
 
402
  } // finished unpacking
 
403
 
 
404
#ifndef NDEBUG
 
405
  len45 = element_list.size();
 
406
  for(unsigned i=0; i<len45; i++)
 
407
    for(unsigned j=0; j<len45; j++)
 
408
      if(i!=j)
 
409
        assert(element_list[i] != element_list[j]);
 
410
#endif
 
411
  
 
412
  ECHO("Halo2 created.");
 
413
  return; 
 
414
 
 
415
}
 
416
  
 
417
 
 
418
 
 
419
 
 
420
 
 
421
 
 
422
 
 
423
 
 
424
 
 
425
 
 
426
 
 
427