← Back to branch summary

~fluidity-core/fluidity/embedded_models

« back to all changes in this revision

Viewing changes to libadaptivity/load_balance/src/fluidity_sam.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/fluidity_sam.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
//
 
29
// Explination of sam options
 
30
//
 
31
 
 
32
// options[0]:  m // Target number of partitions. If this is 0, then
 
33
//                // the number of ranks in MPI_COMM_WORLD is assumed.
 
34
// options[1]:  1 // clean repartitioning into nparts
 
35
//              2 // re-partitioning (diffusive) - local optimization
 
36
//              3 // re-partitioning (directed diffusion) - global optimization
 
37
//              4 // re-partitioning (repartition-remapping)
 
38
// options[2]:  1 // homogenious processing power
 
39
//              2 // hetrogenious processing power
 
40
// options[3]:  1 // no node weights
 
41
//              2 // calculate edge weights based on the projected 
 
42
//                   density of nodes in the future
 
43
// options[4]:  1 // no edge weights
 
44
//              2 // calculate edge weights optimizing for adaptivity
 
45
//              3 // calculate edge weights based on curvature of fields
 
46
// options[5]:  1 // simple mesh
 
47
//              2 // mixed formulation
 
48
#include "confdefs.h"
 
49
 
 
50
#include <map>
 
51
#include <string>
 
52
#include <vector>
 
53
 
 
54
#include "c++debug.h"
 
55
#include "Mesh.h"
 
56
#include "samtypes.h"
 
57
 
 
58
using namespace std;
 
59
 
 
60
#define sam_fc F77_FUNC(sam, SAM)
 
61
#define sam_init_fc F77_FUNC(sam_init_c, SAM_INIT_C)
 
62
#define sam_migrate_fc F77_FUNC(sam_migrate_c, SAM_MIGRATE_C)
 
63
#define sam_query_fc F77_FUNC(sam_query_c, SAM_QUERY_C)
 
64
#define sam_export_mesh_fc F77_FUNC(sam_export_mesh_c, SAM_EXPORT_MESH_C)
 
65
#define sam_export_halo_fc F77_FUNC(sam_export_halo_c, SAM_EXPORT_HALO_C)
 
66
#define sam_export_phalo_fc F77_FUNC(sam_export_phalo_c, SAM_EXPORT_PHALO_C)
 
67
#define sam_cleanup_fc F77_FUNC(sam_cleanup_c, SAM_CLEANUP_C)
 
68
#define sam_add_field_fc F77_FUNC(sam_add_field_c, SAM_ADD_FIELD_C)
 
69
#define sam_pop_field_fc F77_FUNC(sam_pop_field_c, SAM_POP_FIELD_C)
 
70
#define sam_export_node_ownership_fc F77_FUNC(sam_export_node_ownership_c, SAM_EXPORT_NODE_OWNERSHIP_C)
 
71
 
 
72
// "Mesh" is the principle object that does everything for SAM.
 
73
Mesh *mesh = NULL;
 
74
vector<int> decomp_opts;
 
75
vector<int> noddom;
 
76
  
 
77
extern "C"{  
 
78
  void sam_init_fc(
 
79
                   const int* dim, const int *NNodes,  const int *NElems, const int *NSElems, 
 
80
                   const int GATHER[], const int ATOSEN[],
 
81
                   const int SCATER[], const int ATOREC[],
 
82
                   const int* ncolga, const int* nscate, const int* nprocs,
 
83
                   const int ENLIST[], const int *nloc,
 
84
                   const int SNLIST[], const int SURFID[], const int *snloc,
 
85
                   const samfloat_t NODX[], const samfloat_t NODY[], const samfloat_t NODZ[],
 
86
                   samfloat_t Metric[],     const samfloat_t FIELDS[], const int *NFIELDS, 
 
87
                   const int options[], const samfloat_t *fxnl_tol){
 
88
#ifdef HAVE_MPI
 
89
    if(mesh != NULL)
 
90
    {
 
91
      delete mesh;
 
92
      mesh = NULL;
 
93
    }
 
94
    mesh = new Mesh;
 
95
 
 
96
    CHECK( options[5] );
 
97
    mesh->mixed_formulation(options[5] == 2);
 
98
  
 
99
  
 
100
    if( mesh->mixed_formulation() )
 
101
      ECHO("Mixed-formulation capabilities have been enabled.");
 
102
  
 
103
    // Slurp in Fluiditys' mesh
 
104
    ECHO("Importing mesh");
 
105
    mesh->import_fluidity(*dim, *NNodes, *NElems, *NSElems, *NNodes,
 
106
                         *NFIELDS, NODX, NODY, NODZ, Metric, FIELDS,
 
107
                         ENLIST, *nloc,
 
108
                         SNLIST, SURFID, *snloc,
 
109
                         ATOREC, SCATER,
 
110
                         ATOSEN, GATHER);
 
111
    ECHO("Mesh imported");
 
112
 
 
113
    // Build options.
 
114
    decomp_opts.resize(5);
 
115
    for(int i=0;i<5;i++)
 
116
      decomp_opts[i] = options[i];
 
117
 
 
118
    mesh->set_functional_tol(*fxnl_tol);
 
119
  
 
120
    if( (options[2] == 2)&&(options[3] == 2) ){
 
121
      cout<<"WARNING: It has been requested to load balance the mesh on a system with "          
 
122
          << "heterogeneous point-to-point bandwidth and heterogeneous processing power."
 
123
          << "This is a non-trivial request. SAM promises to do his best, but don't "     
 
124
          << "expect wonders."<< endl;
 
125
    }
 
126
#endif
 
127
    return;
 
128
  }
 
129
 
 
130
  void sam_migrate_fc(void)
 
131
  {
 
132
#ifdef HAVE_MPI
 
133
    assert(mesh != NULL);
 
134
    
 
135
    // Graph partitioning.
 
136
    ECHO("Partitioning mesh...");
 
137
    noddom = mesh->decomp( decomp_opts );
 
138
    ECHO("...partitioned.");
 
139
  
 
140
    // Migrate Mesh according to noddom.
 
141
    ECHO("Migrating mesh...");
 
142
    mesh->migrate( noddom );
 
143
    ECHO("...migrated.");
 
144
  
 
145
    if(mesh->mixed_formulation()){
 
146
      mesh->invent_pressure_mesh();
 
147
      mesh->formHalo2();
 
148
    }
 
149
#endif
 
150
    return;
 
151
  }
 
152
 
 
153
  void sam_query_fc(int* NONODS, int* TOTELE, int* STOTEL,
 
154
                    int* ncolga, int* nscate, int* pncolga, int* pnscate)
 
155
  {
 
156
#ifdef HAVE_MPI
 
157
    assert(mesh != NULL);
 
158
 
 
159
    mesh->do_element_headcount();
 
160
    *NONODS = mesh->node_list.size();
 
161
    *TOTELE = mesh->num_elements("volume");
 
162
    *STOTEL = mesh->num_elements("surface");
 
163
 
 
164
    *ncolga = mesh->get_ncolga();
 
165
    *nscate = mesh->get_nscate();
 
166
 
 
167
    if(mesh->mixed_formulation())
 
168
    {
 
169
      *pncolga = mesh->get_pncolga();
 
170
      *pnscate = mesh->get_pnscate();
 
171
    }
 
172
    else
 
173
    {
 
174
      *pncolga = -1;
 
175
      *pnscate = -1;
 
176
    }
 
177
#endif
 
178
    return;
 
179
  }
 
180
 
 
181
  void sam_export_mesh_fc(int* nonods, int* totele, int* stotel, int* nloc, int* snloc,
 
182
      samfloat_t* NODX, samfloat_t* NODY, samfloat_t* NODZ, int* ENLIST, int* SENLIST, int* SURFID)
 
183
  {
 
184
#ifdef HAVE_MPI
 
185
    assert(mesh != NULL);
 
186
 
 
187
    // Positions field
 
188
    {
 
189
      int i = 0;
 
190
      for(deque<Node>::iterator in=mesh->node_list.begin(); in != mesh->node_list.end(); ++in){
 
191
        switch(in->get_size_x()){
 
192
          case 3:
 
193
            NODZ[i] = in->get_z();
 
194
          case 2:
 
195
            NODY[i] = in->get_y();
 
196
          case 1:
 
197
            NODX[i] = in->get_x();
 
198
            break;
 
199
          default:
 
200
            ERROR("Invalid dimension");
 
201
        }
 
202
        i++;
 
203
      }
 
204
    }
 
205
    
 
206
    int num_elements = mesh->element_list.size();
 
207
    vector<int> elm_owner(num_elements);
 
208
    for(int i=0;i<num_elements;i++){
 
209
      vector<unn_t> enl(mesh->element_list[i].get_enlist());
 
210
      
 
211
      vector<unn_t>::iterator it=enl.begin();
 
212
      int gnn = mesh->unn2gnn(*it);
 
213
      elm_owner[i] = mesh->node_list[gnn].get_future_owner();
 
214
      
 
215
      for(;it!=enl.end();++it){
 
216
        gnn = mesh->unn2gnn(*it);
 
217
        elm_owner[i] = min(elm_owner[i], (int)mesh->node_list[gnn].get_future_owner());
 
218
      }
 
219
    }
 
220
 
 
221
    int MyRank = MPI::COMM_WORLD.Get_rank();
 
222
 
 
223
    { // Compress volume element-node lists
 
224
      int i = 0;
 
225
      vector<int> halo_elements;
 
226
      for(int e=0;e<num_elements;e++){
 
227
        unsigned char type = mesh->element_list[e].get_flags();
 
228
        if( type & ELM_VOLUME ){
 
229
          vector<unn_t> enl(mesh->element_list[e].get_enlist());
 
230
                    
 
231
          if(elm_owner[e]==MyRank){
 
232
            for(vector<unn_t>::const_iterator it=enl.begin(); it!=enl.end(); ++it){
 
233
              ENLIST[i++] = mesh->unn2gnn(*it) + 1;
 
234
            }
 
235
          }else{
 
236
            for(vector<unn_t>::const_iterator it=enl.begin(); it!=enl.end(); ++it){
 
237
              halo_elements.push_back(mesh->unn2gnn(*it) + 1);
 
238
            }
 
239
          }
 
240
        }
 
241
      }
 
242
      for(vector<int>::const_iterator it=halo_elements.begin();it!=halo_elements.end();++it){
 
243
        ENLIST[i++]=*it;
 
244
      }
 
245
    }
 
246
 
 
247
    { 
 
248
      // Compress surface element-node list and write the surface id's.    
 
249
      mesh->do_element_headcount();
 
250
#ifndef NDEBUG
 
251
      int NewNSElems = mesh->num_elements("surface");
 
252
#endif
 
253
      int pos=0;
 
254
      int i = 0;
 
255
      vector<int> halo_elements, halo_element_ids;
 
256
      for(int e=0;e<num_elements;e++){
 
257
        unsigned char type = mesh->element_list[e].get_flags();
 
258
        if( type & ELM_SURFACE ){
 
259
          assert(pos<NewNSElems);
 
260
          
 
261
          vector<unn_t> enl( mesh->element_list[e].get_enlist() );
 
262
          const vector<int>& ifields = mesh->element_list[e].get_ifields();
 
263
          if(elm_owner[e]==MyRank){
 
264
            for(vector<unn_t>::const_iterator it = enl.begin(); it != enl.end(); ++it){
 
265
              SENLIST[i++] = mesh->unn2gnn(*it) + 1;
 
266
            }
 
267
            SURFID[pos++] = ifields[0];
 
268
          }else{
 
269
            for(vector<unn_t>::const_iterator it = enl.begin(); it != enl.end(); ++it){
 
270
              halo_elements.push_back(mesh->unn2gnn(*it) + 1);
 
271
            }
 
272
            halo_element_ids.push_back(ifields[0]);            
 
273
          }
 
274
        }
 
275
      }
 
276
      for(vector<int>::const_iterator it=halo_elements.begin();it!=halo_elements.end();++it){
 
277
        SENLIST[i++]=*it;
 
278
      }
 
279
      for(vector<int>::const_iterator it=halo_element_ids.begin();it!=halo_element_ids.end();++it){
 
280
        SURFID[pos++]=*it;
 
281
      }
 
282
    }
 
283
#endif
 
284
    return;
 
285
  }
 
286
  
 
287
  void sam_export_halo_fc(int* colgat, int* atosen, int* scater, int* atorec, const int* ncolga, const int* nscate, const int* nprocs, int* pnodes, int* nnodes){
 
288
#ifdef HAVE_MPI
 
289
    assert(mesh != NULL);
 
290
    
 
291
    mesh->export_halo(colgat, atosen, scater, atorec, ncolga, nscate, nprocs);
 
292
 
 
293
    *pnodes = mesh->node_list.psize();
 
294
    *nnodes = mesh->node_list.size();
 
295
#endif
 
296
    return;
 
297
  }
 
298
 
 
299
  void sam_export_phalo_fc(int* pcolgat, int* patosen, int* pscater, int* patorec, const int* pncolga, const int* pnscate, const int* nprocs, int* ppnodes, int* pnnodes){
 
300
#ifdef HAVE_MPI
 
301
    assert(mesh != NULL);
 
302
    
 
303
    mesh->export_phalo(pcolgat, patosen, pscater, patorec, pncolga, pnscate, nprocs);
 
304
    *ppnodes = mesh->MFnode_list.psize();
 
305
    *pnnodes = mesh->MFnode_list.size();
 
306
#endif
 
307
    return;
 
308
  }
 
309
  
 
310
  void sam_cleanup_fc(void){
 
311
    if(mesh != NULL){
 
312
      delete mesh;
 
313
      mesh = NULL;
 
314
      
 
315
      decomp_opts.clear();
 
316
    }
 
317
    noddom.clear();
 
318
    
 
319
    return;
 
320
  }
 
321
 
 
322
  void sam_add_field_fc(samfloat_t* field_data, int *nnodes){
 
323
#ifdef HAVE_MPI
 
324
    assert(mesh != NULL);
 
325
 
 
326
    assert(mesh->node_list.size() == (size_t)*nnodes);
 
327
    for(size_t i=0;i<mesh->node_list.size();i++){
 
328
      mesh->node_list[i].append_field(field_data + i, 1);
 
329
    }
 
330
#endif
 
331
    return;
 
332
  }
 
333
  
 
334
  void sam_pop_field_fc(samfloat_t* field_data, int *nnodes){
 
335
#ifdef HAVE_MPI
 
336
    assert(mesh != NULL);
 
337
    
 
338
    assert(mesh->node_list.size() == (size_t)*nnodes);
 
339
    for (size_t i=0;i<mesh->node_list.size();i++){
 
340
      field_data[i] = mesh->node_list[i].pop_field();
 
341
    }
 
342
#endif
 
343
    return;
 
344
  }  
 
345
  
 
346
  void sam_export_node_ownership_fc(int* node_ownership, int* nnodes){
 
347
    assert(noddom.size() == *nnodes);
 
348
      
 
349
    for(int i = 0;i < *nnodes;i++){
 
350
      node_ownership[i] = noddom[i];
 
351
    }
 
352
    
 
353
    return;
 
354
  }
 
355
}