1 // Copyright (c) Lawrence Livermore National Security, LLC and other Conduit
2 // Project developers. See top-level LICENSE AND COPYRIGHT files for dates and
3 // other details. No copyright assignment is required to contribute to Conduit.
4
5 //-----------------------------------------------------------------------------
6 ///
7 /// file: conduit_blueprint_mpi_mesh_parmetis.cpp
8 ///
9 //-----------------------------------------------------------------------------
10
11 //-----------------------------------------------------------------------------
12 // conduit includes
13 //-----------------------------------------------------------------------------
14 #include "conduit_blueprint_mesh.hpp"
15 #include "conduit_blueprint_mesh_utils.hpp"
16 #include "conduit_blueprint_mpi_mesh.hpp"
17 #include "conduit_blueprint_mpi_mesh_parmetis.hpp"
18 #include "conduit_blueprint_o2mrelation.hpp"
19 #include "conduit_blueprint_o2mrelation_iterator.hpp"
20
21 #include "conduit_relay_mpi.hpp"
22
23 #include <parmetis.h>
24
25 //-----------------------------------------------------------------------------
26 // -- begin conduit --
27 //-----------------------------------------------------------------------------
28 namespace conduit
29 {
30
31 //-----------------------------------------------------------------------------
32 // -- begin conduit::blueprint --
33 //-----------------------------------------------------------------------------
34 namespace blueprint
35 {
36
37 //-----------------------------------------------------------------------------
38 // -- begin conduit::blueprint::mpi --
39 //-----------------------------------------------------------------------------
40 namespace mpi
41 {
42
43 //-----------------------------------------------------------------------------
44 // -- begin conduit::blueprint::mesh --
45 //-----------------------------------------------------------------------------
46
47 namespace mesh
48 {
49 //-----------------------------------------------------------------------------
50
51 //-----------------------------------------------------------------------------
52
53 //-----------------------------------------------------------------------------
54 //-----------------------------------------------------------------------------
55 //-- Map Parmetis Types (idx_t and real_t) to conduit dtype ids
56 //-----------------------------------------------------------------------------
57 //-----------------------------------------------------------------------------
58
59 // check our assumptions
60 static_assert(IDXTYPEWIDTH != 32 || IDXTYPEWIDTH != 64,
61 "Metis idx_t is not 32 or 64 bits");
62
63 static_assert(REALTYPEWIDTH != 32 || REALTYPEWIDTH != 64,
64 "Metis real_t is not 32 or 64 bits");
65
66 //-----------------------------------------------------------------------------
67 //-----------------------------------------------------------------------------
68 // IDXTYPEWIDTH and REALTYPEWIDTH are metis type defs
69 //-----------------------------------------------------------------------------
70 //-----------------------------------------------------------------------------
71
72 //-----------------------------------------------------------------------------
73 index_t
metis_idx_t_to_conduit_dtype_id()74 metis_idx_t_to_conduit_dtype_id()
75 {
76 #if IDXTYPEWIDTH == 64
77 // 64 bits
78 // int64
79 return conduit::DataType::INT64_ID;
80 #else
81 // 32 bits
82 // int32
83 return conduit::DataType::INT32_ID;
84 #endif
85 }
86
87 //-----------------------------------------------------------------------------
88 index_t
metis_real_t_t_to_conduit_dtype_id()89 metis_real_t_t_to_conduit_dtype_id()
90 {
91 #if REALTYPEWIDTH == 64
92 // 64 bits
93 // float64
94 return conduit::DataType::FLOAT64_ID;
95 #else
96 // 32 bits
97 // float32
98 return conduit::DataType::FLOAT32_ID;
99 #endif
100 }
101
102
103 //-----------------------------------------------------------------------------
104 // NOTE: this is generally useful, it should be added to mpi::mesh
105 //
106 // supported options:
107 // topology: {string}
108 // field_prefix: {string}
generate_global_element_and_vertex_ids(conduit::Node & mesh,const Node & options,MPI_Comm comm)109 void generate_global_element_and_vertex_ids(conduit::Node &mesh,
110 const Node &options,
111 MPI_Comm comm)
112 {
113 // TODO: Check of dest fields already exist, if they do error
114
115
116 int par_rank = conduit::relay::mpi::rank(comm);
117 int par_size = conduit::relay::mpi::size(comm);
118
119 index_t local_num_doms = ::conduit::blueprint::mesh::number_of_domains(mesh);
120 index_t global_num_doms = number_of_domains(mesh,comm);
121
122 if(global_num_doms == 0)
123 {
124 return;
125 }
126
127 std::vector<Node*> domains;
128 ::conduit::blueprint::mesh::domains(mesh,domains);
129
130 // parse options
131 std::string topo_name = "";
132 std::string field_prefix = "";
133 if( options.has_child("topology") )
134 {
135 topo_name = options["topology"].as_string();
136 }
137 else
138 {
139 // TOOD: IMP find the first topo name on a rank with data
140 // for now, just grab first topo
141 const Node &dom_topo = domains[0]->fetch("topologies")[0];
142 topo_name = dom_topo.name();
143 }
144
145 if( options.has_child("field_prefix") )
146 {
147 field_prefix = options["field_prefix"].as_string();
148 }
149
150 // count all local elements + verts and create offsets
151 uint64 local_total_num_eles=0;
152 uint64 local_total_num_verts=0;
153
154 Node local_info;
155 // per domain local info books
156 local_info["num_verts"].set(DataType::uint64(local_num_doms));
157 local_info["num_eles"].set(DataType::uint64(local_num_doms));
158 local_info["verts_offsets"].set(DataType::uint64(local_num_doms));
159 local_info["eles_offsets"].set(DataType::uint64(local_num_doms));
160
161 uint64_array local_num_verts = local_info["num_verts"].value();
162 uint64_array local_num_eles = local_info["num_eles"].value();
163
164 uint64_array local_vert_offsets = local_info["verts_offsets"].value();
165 uint64_array local_ele_offsets = local_info["eles_offsets"].value();
166
167
168 for(size_t local_dom_idx=0; local_dom_idx < domains.size(); local_dom_idx++)
169 {
170 Node &dom = *domains[local_dom_idx];
171 // we do need to make sure we have the requested topo
172 if(dom["topologies"].has_child(topo_name))
173 {
174 // get the topo node
175 const Node &dom_topo = dom["topologies"][topo_name];
176 // get the number of elements in the topo
177 local_num_eles[local_dom_idx] = blueprint::mesh::utils::topology::length(dom_topo);
178 local_ele_offsets[local_dom_idx] = local_total_num_eles;
179 local_total_num_eles += local_num_eles[local_dom_idx];
180
181 // get the coordset that the topo refs
182 const Node &dom_cset = dom["coordsets"][dom_topo["coordset"].as_string()];
183 // get the number of points in the coordset
184 // THIS RETURNS ZERO:
185 //local_num_verts[local_dom_idx] = blueprint::mesh::utils::coordset::length(dom_cset);
186 // so we are using this:
187 local_num_verts[local_dom_idx] = dom_cset["values/x"].dtype().number_of_elements();
188 local_vert_offsets[local_dom_idx] = local_total_num_verts;
189 local_total_num_verts += local_num_verts[local_dom_idx];
190 }
191 }
192
193 // calc per MPI task offsets using
194 // local_total_num_verts
195 // local_total_num_eles
196
197 // first count verts
198 Node max_local, max_global;
199 max_local.set(DataType::uint64(par_size));
200 max_global.set(DataType::uint64(par_size));
201
202 uint64_array max_local_vals = max_local.value();
203 uint64_array max_global_vals = max_global.value();
204
205 max_local_vals[par_rank] = local_total_num_verts;
206
207 relay::mpi::max_all_reduce(max_local, max_global, comm);
208
209
210 index_t global_verts_offset = 0;
211 for(index_t i=0; i< par_rank; i++ )
212 {
213 global_verts_offset += max_global_vals[i];
214 }
215
216 // reset our buffers
217 for(index_t i=0; i< par_size; i++ )
218 {
219 max_local_vals[i] = 0;
220 max_global_vals[i] = 0;
221 }
222
223 // now count eles
224 max_local_vals[par_rank] = local_total_num_eles;
225
226 relay::mpi::max_all_reduce(max_local, max_global, comm);
227
228 index_t global_eles_offset = 0;
229 for(index_t i=0; i< par_rank; i++ )
230 {
231 global_eles_offset += max_global_vals[i];
232 }
233
234 // we now have our offsets, we can create output fields on each local domain
235 for(size_t local_dom_idx=0; local_dom_idx < domains.size(); local_dom_idx++)
236 {
237 Node &dom = *domains[local_dom_idx];
238 // we do need to make sure we have the requested topo
239 if(dom["topologies"].has_child(topo_name))
240 {
241 Node &verts_field = dom["fields"][field_prefix + "global_vertex_ids"];
242 verts_field["association"] = "vertex";
243 verts_field["topology"] = topo_name;
244 verts_field["values"].set(DataType::int64(local_num_verts[local_dom_idx]));
245
246 int64 vert_base_idx = global_verts_offset + local_vert_offsets[local_dom_idx];
247
248 int64_array vert_ids_vals = verts_field["values"].value();
249 for(uint64 i=0; i < local_num_verts[local_dom_idx]; i++)
250 {
251 vert_ids_vals[i] = i + vert_base_idx;
252 }
253
254 // NOTE: VISIT BP DOESNT SUPPORT UINT64!!!!
255 Node &eles_field = dom["fields"][field_prefix + "global_element_ids"];
256 eles_field["association"] = "element";
257 eles_field["topology"] = topo_name;
258 eles_field["values"].set(DataType::int64(local_num_eles[local_dom_idx]));
259
260 int64 ele_base_idx = global_eles_offset + local_ele_offsets[local_dom_idx];
261
262 int64_array ele_ids_vals = eles_field["values"].value();
263 for(uint64 i=0; i < local_num_eles[local_dom_idx]; i++)
264 {
265 ele_ids_vals[i] = i + ele_base_idx;
266 }
267 }
268 }
269 }
270
271 //-----------------------------------------------------------------------------
generate_partition_field(conduit::Node & mesh,MPI_Comm comm)272 void generate_partition_field(conduit::Node &mesh,
273 MPI_Comm comm)
274 {
275 Node opts;
276 generate_partition_field(mesh,opts,comm);
277 }
278
279 //-----------------------------------------------------------------------------
generate_partition_field(conduit::Node & mesh,const conduit::Node & options,MPI_Comm comm)280 void generate_partition_field(conduit::Node &mesh,
281 const conduit::Node &options,
282 MPI_Comm comm)
283 {
284 generate_global_element_and_vertex_ids(mesh,
285 options,
286 comm);
287
288 int par_rank = conduit::relay::mpi::rank(comm);
289 int par_size = conduit::relay::mpi::size(comm);
290
291 index_t global_num_doms = number_of_domains(mesh,comm);
292
293 if(global_num_doms == 0)
294 {
295 return;
296 }
297
298 std::vector<Node*> domains;
299 ::conduit::blueprint::mesh::domains(mesh,domains);
300
301 // parse options
302 std::string topo_name = "";
303 std::string field_prefix = "";
304 idx_t nparts = (idx_t)global_num_doms;
305
306 if( options.has_child("topology") )
307 {
308 topo_name = options["topology"].as_string();
309 }
310 else
311 {
312 // TOOD: IMP find the first topo name on a rank with data
313 // for now, just grab first topo
314 const Node &dom_topo = domains[0]->fetch("topologies")[0];
315 topo_name = dom_topo.name();
316 }
317 idx_t ncommonnodes;
318 if ( options.has_child("parmetis_ncommonnodes") )
319 {
320 ncommonnodes = options["parmetis_ncommonnodes"].as_int();
321 }
322 else
323 {
324 // in 2D, zones adjacent if they share 2 nodes (edge)
325 // in 3D, zones adjacent if they share 3 nodes (plane)
326 std::string coordset_name
327 = domains[0]->fetch(std::string("topologies/") + topo_name + "/coordset").as_string();
328 const Node& coordset = domains[0]->fetch(std::string("coordsets/") + coordset_name);
329 ncommonnodes = conduit::blueprint::mesh::coordset::dims(coordset);
330 }
331
332 if( options.has_child("field_prefix") )
333 {
334 field_prefix = options["field_prefix"].as_string();
335 }
336
337 if( options.has_child("partitions") )
338 {
339 nparts = (idx_t) options["partitions"].to_int64();
340 }
341 // TODO: Should this be an error or use default (discuss more)?
342 // else
343 // {
344 // CONDUIT_ERROR("Missing required option in generate_partition_field(): "
345 // << "expected \"partitions\" field with number of partitions.");
346 // }
347
348 // we now have global element and vertex ids
349 // we just need to do some counting and then
350 // traverse our topo to convert this info to parmetis input
351
352 // we need the total number of local eles
353 // the total number of element to vers entries
354
355 index_t local_total_num_eles =0;
356 index_t local_total_ele_to_verts_size = 0;
357
358 for(size_t local_dom_idx=0; local_dom_idx < domains.size(); local_dom_idx++)
359 {
360 Node &dom = *domains[local_dom_idx];
361 // we do need to make sure we have the requested topo
362 if(dom["topologies"].has_child(topo_name))
363 {
364 // get the topo node
365 const Node &dom_topo = dom["topologies"][topo_name];
366 // get the number of elements in the topo
367 local_total_num_eles += blueprint::mesh::utils::topology::length(dom_topo);
368
369 Node topo_offsets;
370 blueprint::mesh::topology::unstructured::generate_offsets(dom_topo, topo_offsets);
371
372 // TODO:
373 // for unstructured we need to do shape math, for unif/rect/struct
374 // we need to do implicit math
375
376 // for unstructured poly:
377 // add up all the sizes, don't use offsets?
378 uint64_accessor sizes_vals = dom_topo["elements/sizes"].as_uint64_accessor();
379 for(index_t i=0; i < sizes_vals.number_of_elements();i++)
380 {
381 local_total_ele_to_verts_size += sizes_vals[i];
382 }
383
384 }
385 }
386
387 // reminder:
388 // idx_t eldist[] = {0, 3, 4};
389 //
390 // idx_t eptr[] = {0,4,8,12};
391 //
392 // idx_t eind[] = {0,1,3,4,
393 // 1,2,4,5,
394 // 3,4,6,7};
395
396 Node parmetis_params;
397 // eldist tells how many elements there are per mpi task,
398 // it will be size par_size + 1
399 parmetis_params["eldist"].set(DataType(metis_idx_t_to_conduit_dtype_id(),
400 par_size+1));
401 // eptr holds the offsets to the start of each element's
402 // vertex list
403 // size == total number of local elements (we counted this above)
404 parmetis_params["eptr"].set(DataType(metis_idx_t_to_conduit_dtype_id(),
405 local_total_num_eles+1));
406 // eind holds, for each element, a list of vertex ids
407 // (we also counted this above)
408 parmetis_params["eind"].set(DataType(metis_idx_t_to_conduit_dtype_id(),
409 local_total_ele_to_verts_size));
410
411 // output array, size of local num elements
412 parmetis_params["part"].set(DataType(metis_idx_t_to_conduit_dtype_id(),
413 local_total_num_eles));
414
415
416 // first lets get eldist setup:
417 // eldist[0] = 0,\
418 // eldist[1] == # of elements on rank 0-
419 // eldist[2] == # of elemens on rank 0 + rank 1
420 // ...
421 // eldist[n] == # of total elements
422 //
423 Node el_counts;
424 el_counts["local"] = DataType(metis_idx_t_to_conduit_dtype_id(),
425 par_size);
426 el_counts["global"] = DataType(metis_idx_t_to_conduit_dtype_id(),
427 par_size);
428
429 idx_t *el_counts_local_vals = el_counts["local"].value();
430 idx_t *el_counts_global_vals = el_counts["global"].value();
431 el_counts_local_vals[par_rank] = local_total_num_eles;
432 relay::mpi::max_all_reduce(el_counts["local"], el_counts["global"], comm);
433
434 // prefix sum to set eldist
435 idx_t *eldist_vals = parmetis_params["eldist"].value();
436 eldist_vals[0] = 0;
437 for(size_t i=0;i<(size_t)par_size;i++)
438 {
439 eldist_vals[i+1] = eldist_vals[i] + el_counts_global_vals[i];
440 }
441
442 idx_t *eptr_vals = parmetis_params["eptr"].value();
443 idx_t *eind_vals = parmetis_params["eind"].value();
444
445 // now elptr == prefix sum of the sizes
446 // (note: the offsets don't matter for elptr b/c we are creating a compact
447 // rep for parmetis)
448 //
449 // and eind == look up of global vertex id
450 size_t eptr_idx=0;
451 size_t eind_idx=0;
452 idx_t curr_offset = 0;
453 for(size_t local_dom_idx=0; local_dom_idx < domains.size(); local_dom_idx++)
454 {
455 Node &dom = *domains[local_dom_idx];
456 // we do need to make sure we have the requested topo
457 if(dom["topologies"].has_child(topo_name))
458 {
459 // get the topo node
460 Node &dom_topo = dom["topologies"][topo_name];
461 const Node &dom_g_vert_ids = dom["fields"][field_prefix + "global_vertex_ids"]["values"];
462
463 // for unstruct poly: use sizes
464 uint64_accessor sizes_vals = dom_topo["elements/sizes"].as_uint64_accessor();
465 for(index_t i=0; i < sizes_vals.number_of_elements(); i++)
466 {
467 eptr_vals[eptr_idx] = curr_offset;
468 curr_offset += sizes_vals[i];
469 eptr_idx++;
470 }
471 // add last offset
472 eptr_vals[eptr_idx] = curr_offset;
473
474 int64_accessor global_vert_ids = dom_g_vert_ids.as_int64_accessor();
475 // for each element:
476 // loop over each local vertex, and use global vert map to add and entry to eind
477
478 uint64_accessor conn_vals = dom_topo["elements/connectivity"].as_uint64_accessor();
479
480 o2mrelation::O2MIterator o2miter(dom_topo["elements"]);
481 while(o2miter.has_next(conduit::blueprint::o2mrelation::ONE))
482 {
483 o2miter.next(conduit::blueprint::o2mrelation::ONE);
484 o2miter.to_front(conduit::blueprint::o2mrelation::MANY);
485 while(o2miter.has_next(conduit::blueprint::o2mrelation::MANY))
486 {
487 o2miter.next(conduit::blueprint::o2mrelation::MANY);
488 const index_t local_vert_id = o2miter.index(conduit::blueprint::o2mrelation::DATA);
489 // get the conn
490 eind_vals[eind_idx] = (idx_t) global_vert_ids[conn_vals[local_vert_id]];
491 eind_idx++;
492 }
493 }
494 }
495 }
496
497 idx_t wgtflag = 0; // weights are NULL
498 idx_t numflag = 0; // C-style numbering
499 idx_t ncon = 1; // the number of weights per vertex
500 // equal weights for each proc
501 std::vector<real_t> tpwgts(nparts, 1.0/nparts);
502 real_t ubvec = 1.050000;
503
504 // options == extra output
505 idx_t parmetis_opts[] = {1,
506 PARMETIS_DBGLVL_TIME |
507 PARMETIS_DBGLVL_INFO |
508 PARMETIS_DBGLVL_PROGRESS |
509 PARMETIS_DBGLVL_REFINEINFO |
510 PARMETIS_DBGLVL_MATCHINFO |
511 PARMETIS_DBGLVL_RMOVEINFO |
512 PARMETIS_DBGLVL_REMAP,
513 0};
514 // outputs
515 idx_t edgecut = 0; // will hold # of cut edges
516
517 // output array, size of local num elements
518 parmetis_params["part"].set(DataType(metis_idx_t_to_conduit_dtype_id(),
519 local_total_num_eles));
520 idx_t *part_vals = parmetis_params["part"].value();
521
522 int parmetis_res = ParMETIS_V3_PartMeshKway(eldist_vals,
523 eptr_vals,
524 eind_vals,
525 NULL,
526 &wgtflag,
527 &numflag,
528 &ncon,
529 &ncommonnodes,
530 &nparts,
531 tpwgts.data(),
532 &ubvec,
533 parmetis_opts,
534 &edgecut,
535 part_vals,
536 &comm);
537
538 index_t part_vals_idx=0;
539 // create output field with part result
540 for(size_t local_dom_idx=0; local_dom_idx < domains.size(); local_dom_idx++)
541 {
542 Node &dom = *domains[local_dom_idx];
543 // we do need to make sure we have the requested topo
544 if(dom["topologies"].has_child(topo_name))
545 {
546 // get the topo node
547 const Node &dom_topo = dom["topologies"][topo_name];
548 // get the number of elements in the topo
549 index_t dom_num_eles = blueprint::mesh::utils::topology::length(dom_topo);
550 // for unstrcut we need to do shape math, for unif/rect/struct
551 // we need to do implicit math
552
553 // NOTE: VISIT BP DOESNT SUPPORT UINT64!!!!
554 Node &part_field = dom["fields"][field_prefix + "parmetis_result"];
555 part_field["association"] = "element";
556 part_field["topology"] = topo_name;
557 part_field["values"].set(DataType::int64(dom_num_eles));
558
559 int64_array part_field_vals = part_field["values"].value();
560 for(index_t i=0; i < dom_num_eles; i++)
561 {
562 part_field_vals[i] = part_vals[part_vals_idx];
563 part_vals_idx++;
564 }
565 }
566 }
567
568 }
569
570
571 //-----------------------------------------------------------------------------
572 }
573 //-----------------------------------------------------------------------------
574 // -- end conduit::blueprint::mpi::mesh --
575 //-----------------------------------------------------------------------------
576
577
578 }
579 //-----------------------------------------------------------------------------
580 // -- end conduit::blueprint::mpi --
581 //-----------------------------------------------------------------------------
582
583
584 }
585 //-----------------------------------------------------------------------------
586 // -- end conduit::blueprint --
587 //-----------------------------------------------------------------------------
588
589 }
590 //-----------------------------------------------------------------------------
591 // -- end conduit:: --
592 //-----------------------------------------------------------------------------
593
594