1 /*============================================================================
2 * Handle boxes aligned with Cartesian axes.
3 *===========================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *---------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <float.h>
35 #include <limits.h>
36 #include <math.h>
37 #include <stdlib.h>
38 #include <string.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *---------------------------------------------------------------------------*/
43
44 #include "bft_mem.h"
45 #include "bft_printf.h"
46
47 /*----------------------------------------------------------------------------
48 * Header for the current file
49 *---------------------------------------------------------------------------*/
50
51 #include "fvm_box.h"
52 #include "fvm_box_priv.h"
53
54 #include "cs_all_to_all.h"
55
56 /*---------------------------------------------------------------------------*/
57
58 BEGIN_C_DECLS
59
60 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
61
62 /*=============================================================================
63 * Local Macro definitions
64 *===========================================================================*/
65
66 /*============================================================================
67 * Type and structure definitions
68 *===========================================================================*/
69
70 /*============================================================================
71 * Private function definitions
72 *===========================================================================*/
73
74 #if defined(HAVE_MPI)
75
76 /*----------------------------------------------------------------------------
77 * Display a histogram on leaves associated to the boxes and several
78 * other pieces of information (min, max, ...)
79 *
80 * parameters:
81 * distrib <-- pointer to the fvm_box_distrib_t structure
82 * n_quantiles <-> number of quantiles requested (or NULL);
83 * may be lower upon return
84 * quantile_start --> start of quantile (size: n_quantiles + 1)
85 * n_quantile_boxes --> number of boxes per quantile (size: n_quantiles)
86 * imbalance --> distribution imbalance measure (or NULL)
87 * n_ranks --> number of ranks with boxes (or NULL)
88 * comm <-- associated MPI communicator
89 *---------------------------------------------------------------------------*/
90
91 static void
_get_distrib_statistics(const fvm_box_distrib_t * distrib,cs_lnum_t * n_quantiles,cs_lnum_t quantile_start[],cs_lnum_t n_quantile_boxes[],double * imbalance,int * n_ranks,MPI_Comm comm)92 _get_distrib_statistics(const fvm_box_distrib_t *distrib,
93 cs_lnum_t *n_quantiles,
94 cs_lnum_t quantile_start[],
95 cs_lnum_t n_quantile_boxes[],
96 double *imbalance,
97 int *n_ranks,
98 MPI_Comm comm)
99 {
100 cs_lnum_t i, j, k, step, delta, _n_rank_boxes;
101
102 int _n_ranks = 0;
103 cs_lnum_t _min = INT_MAX, _max = 0, gmin = 0, gmax = 0;
104
105 /* Sanity checks */
106
107 assert(distrib != NULL);
108 assert(distrib->index != NULL);
109
110 if (n_quantiles != NULL) {
111
112 cs_lnum_t _n_quantiles = 1;
113
114 /* Get global min and max number of boxes */
115
116 for (i = 0; i < distrib->n_ranks; i++) {
117
118 _n_rank_boxes = distrib->index[i+1] - distrib->index[i];
119 _min = CS_MIN(_min, _n_rank_boxes);
120 _max = CS_MAX(_max, _n_rank_boxes);
121
122 if (_n_rank_boxes > 0)
123 _n_ranks += 1;
124
125 }
126
127 gmin = _min;
128 gmax = _max;
129
130 MPI_Allreduce(&_min, &gmin, 1, CS_MPI_LNUM, MPI_MIN, comm);
131 MPI_Allreduce(&_max, &gmax, 1, CS_MPI_LNUM, MPI_MAX, comm);
132
133 /* Build a histogram for the distribution of boxes */
134
135 delta = gmax - gmin;
136 if (delta < _n_quantiles)
137 _n_quantiles = delta;
138
139 /* Define quantiles */
140
141 step = delta / _n_quantiles;
142 if (delta % _n_quantiles > 0)
143 step++;
144
145 for (i = 0; i < _n_quantiles; i++)
146 quantile_start[i] = gmin + i*step;
147
148 quantile_start[_n_quantiles] = gmax + 1;
149
150 /* Count for quantiles */
151
152 for (j = 0; j < _n_quantiles; j++)
153 n_quantile_boxes[j] = 0;
154
155 if (delta > 0) { /* Loop on boxes */
156
157 for (i = 0; i < distrib->n_ranks; i++) {
158
159 _n_rank_boxes = distrib->index[i+1] - distrib->index[i];
160
161 for (k = 1; k < _n_quantiles; k++) {
162 if (_n_rank_boxes < gmin + k*step)
163 break;
164 }
165 n_quantile_boxes[k-1] += 1;
166
167 } /* End of loop on boxes */
168
169 }
170
171 *n_quantiles = _n_quantiles;
172 }
173
174 /* Set other return values */
175
176 if (imbalance != NULL)
177 *imbalance = distrib->fit;
178
179 if (n_ranks != NULL)
180 *n_ranks = _n_ranks;
181 }
182
183 #endif /* defined(HAVE_MPI) */
184
185 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
186
187 /*============================================================================
188 * Public function definitions
189 *===========================================================================*/
190
191 /*----------------------------------------------------------------------------
192 * Create a set of boxes and initialize it.
193 *
194 * parameters:
195 * dim <-- spatial dimension
196 * normalize <-- 1 if boxes are to be normalized, 0 otherwize
197 * allow_projection <-- if 1, project to lower dimension if all boxes
198 * are cut by the median plane of the set.
199 * n_boxes <-- number of elements to create
200 * box_gnum <-- global numbering of boxes
201 * extents <-- coordinate extents (size: n_boxes*dim*2, as
202 * xmin1, ymin1, .. xmax1, ymax1, ..., xmin2, ...)
203 * comm <-- associated MPI communicator
204 *
205 * returns:
206 * a new allocated pointer to a fvm_box_set_t structure.
207 *---------------------------------------------------------------------------*/
208
209 #if defined(HAVE_MPI)
210 fvm_box_set_t *
fvm_box_set_create(int dim,int normalize,int allow_projection,cs_lnum_t n_boxes,const cs_gnum_t * box_gnum,const cs_coord_t * box_extents,MPI_Comm comm)211 fvm_box_set_create(int dim,
212 int normalize,
213 int allow_projection,
214 cs_lnum_t n_boxes,
215 const cs_gnum_t *box_gnum,
216 const cs_coord_t *box_extents,
217 MPI_Comm comm)
218 #else
219 fvm_box_set_t *
220 fvm_box_set_create(int dim,
221 int normalize,
222 int allow_projection,
223 cs_lnum_t n_boxes,
224 const cs_gnum_t *box_gnum,
225 const cs_coord_t *box_extents)
226 #endif
227 {
228 int j, k;
229 cs_lnum_t i;
230 cs_gnum_t n_g_boxes = n_boxes;
231 cs_coord_t g_min[3], g_max[3], g_extents[6];
232
233 fvm_box_set_t *boxes = NULL;
234
235 /* Get global min/max coordinates */
236
237 #if defined(HAVE_MPI)
238 fvm_morton_get_global_extents(dim, n_boxes, box_extents, g_extents, comm);
239 #else
240 fvm_morton_get_global_extents(dim, n_boxes, box_extents, g_extents);
241 #endif
242
243 for (j = 0; j < 3; j++) {
244 g_min[j] = g_extents[j];
245 g_max[j] = g_extents[j+dim];
246 }
247
248 #if defined(HAVE_MPI)
249
250 if (comm != MPI_COMM_NULL) {
251
252 cs_gnum_t box_max = 0;
253
254 for (i = 0; i < n_boxes; i++)
255 box_max = CS_MAX(box_max, box_gnum[i]);
256
257 MPI_Allreduce(&box_max, &n_g_boxes, 1, CS_MPI_GNUM, MPI_MAX, comm);
258
259 }
260
261 #endif
262
263 /* Allocate box set structure and initialize it */
264
265 BFT_MALLOC(boxes, 1, fvm_box_set_t);
266
267 boxes->dim = dim;
268 boxes->n_boxes = n_boxes;
269 boxes->n_g_boxes = n_g_boxes;
270
271 for (j = 0; j < 3; j++) {
272 boxes->dimensions[j] = j;
273 boxes->gmin[j] = g_min[j];
274 boxes->gmax[j] = g_max[j];
275 }
276
277 boxes->g_num = NULL;
278 boxes->extents = NULL;
279
280 #if defined(HAVE_MPI)
281 boxes->comm = comm;
282 #endif
283
284 /* Optionally allow and detect a layout of lower
285 dimension than the spatial dimension */
286
287 if (allow_projection) {
288
289 double g_mid[3];
290 int proj[] = {1, 1, 1};
291
292 for (j = 0; j < dim; j++)
293 g_mid[j] = (g_min[j] + g_max[j]) * 0.5;
294
295 for (i = 0; i < n_boxes; i++) {
296 for (j = 0; j < dim; j++) {
297 if ( box_extents[i*dim*2 + j] > g_mid[j]
298 || box_extents[i*dim*2 + j+dim] < g_mid[j])
299 proj[j] = 0;
300 }
301 }
302
303 #if defined(HAVE_MPI)
304 if (comm != MPI_COMM_NULL) {
305 int l_proj[3];
306 for (j = 0; j < dim; j++)
307 l_proj[j] = proj[j];
308 MPI_Allreduce(l_proj, proj, dim, MPI_INT, MPI_MIN, comm);
309 }
310 #endif
311
312 boxes->dim = 0;
313 for (j = 0; j < dim; j++) {
314 if (proj[j] == 0) {
315 boxes->dimensions[boxes->dim] = j;
316 boxes->dim += 1;
317 }
318 }
319
320
321 }
322
323 for (j = boxes->dim; j < 3; j++) /* ensure all is initialized */
324 boxes->dimensions[j] = -1;
325
326 /* Now assign values */
327
328 BFT_MALLOC(boxes->g_num, n_boxes, cs_gnum_t);
329 BFT_MALLOC(boxes->extents, n_boxes*boxes->dim*2, cs_coord_t);
330
331 for (i = 0; i < n_boxes; i++) {
332
333 cs_coord_t *_min = boxes->extents + (boxes->dim*2*i);
334 cs_coord_t *_max = _min + boxes->dim;
335
336 boxes->g_num[i] = box_gnum[i];
337
338 for (j = 0; j < boxes->dim; j++) {
339 k = boxes->dimensions[j];
340 _min[j] = box_extents[i*dim*2 + k];
341 _max[j] = box_extents[i*dim*2 + k+dim];
342 assert(_min[j] <= _max[j]);
343 }
344 }
345
346 /* Define the normalized min/max coordinates of the box */
347
348 if (normalize) {
349
350 cs_coord_t d[3], s[3];
351
352 for (j = 0; j < boxes->dim; j++) {
353 k = boxes->dimensions[j];
354 s[j] = g_min[k];
355 d[j] = g_max[k] - g_min[k];
356 }
357
358 for (i = 0; i < n_boxes; i++) {
359
360 cs_coord_t *_min = boxes->extents + (boxes->dim*2*i);
361 cs_coord_t *_max = _min + boxes->dim;
362
363 for (j = 0; j < boxes->dim; j++) {
364 _min[j] = (_min[j] - s[j]) / d[j];
365 _max[j] = (_max[j] - s[j]) / d[j];
366 }
367 }
368
369 }
370
371 /* Return pointer to structure */
372
373 return boxes;
374 }
375
376 /*----------------------------------------------------------------------------
377 * Delete a fvm_box_set_t structure.
378 *
379 * parameters:
380 * boxes <-> pointer to the fvm_box_set_t structure to delete
381 *---------------------------------------------------------------------------*/
382
383 void
fvm_box_set_destroy(fvm_box_set_t ** boxes)384 fvm_box_set_destroy(fvm_box_set_t **boxes)
385 {
386 if (boxes != NULL) {
387
388 fvm_box_set_t *_boxes = *boxes;
389
390 if (_boxes == NULL)
391 return;
392
393 BFT_FREE(_boxes->g_num);
394 BFT_FREE(_boxes->extents);
395 BFT_FREE(_boxes);
396 }
397 }
398
399 /*----------------------------------------------------------------------------
400 * Return the dimension associated with a set of boxes.
401 *
402 * parameters:
403 * boxes <-- pointer to set of boxes
404 *
405 * returns:
406 * associated spatial dimension
407 *---------------------------------------------------------------------------*/
408
409 int
fvm_box_set_get_dim(const fvm_box_set_t * boxes)410 fvm_box_set_get_dim(const fvm_box_set_t *boxes)
411 {
412 int retval = 0;
413
414 if (boxes != NULL)
415 retval = boxes->dim;
416
417 return retval;
418 }
419
420 /*----------------------------------------------------------------------------
421 * Return the local number of boxes in a set.
422 *
423 * parameters:
424 * boxes <-- pointer to set of boxes
425 *
426 * returns:
427 * local number of boxes
428 *---------------------------------------------------------------------------*/
429
430 cs_lnum_t
fvm_box_set_get_size(const fvm_box_set_t * boxes)431 fvm_box_set_get_size(const fvm_box_set_t *boxes)
432 {
433 cs_lnum_t retval = 0;
434
435 if (boxes != NULL)
436 retval = boxes->n_boxes;
437
438 return retval;
439 }
440
441 /*----------------------------------------------------------------------------
442 * Return the global number of boxes in a set.
443 *
444 * parameters:
445 * boxes <-- pointer to set of boxes
446 *
447 * returns:
448 * global number of boxes
449 *---------------------------------------------------------------------------*/
450
451 cs_gnum_t
fvm_box_set_get_global_size(const fvm_box_set_t * boxes)452 fvm_box_set_get_global_size(const fvm_box_set_t *boxes)
453 {
454 cs_gnum_t retval = 0;
455
456 if (boxes != NULL)
457 retval = boxes->n_g_boxes;
458
459 return retval;
460 }
461
462 /*----------------------------------------------------------------------------
463 * Return extents associated with a set of boxes.
464 *
465 * The extents array is organized in the following fashion:
466 * {x_min_0, y_min_0, ..., x_max_0, y_max_0, ...
467 * x_min_n, y_min_n, ..., x_max_n, y_max_n, ...}
468 *
469 * Its size is thus: n_boxes * dim * 2.
470 *
471 * parameters:
472 * boxes <-- pointer to set of boxes
473 *
474 * returns:
475 * pointer to extents array
476 *---------------------------------------------------------------------------*/
477
478 const cs_coord_t *
fvm_box_set_get_extents(fvm_box_set_t * boxes)479 fvm_box_set_get_extents(fvm_box_set_t *boxes)
480 {
481 assert(boxes != NULL);
482
483 return boxes->extents;
484 }
485
486 /*----------------------------------------------------------------------------
487 * Return global numbers associated with a set of boxes.
488 *
489 * parameters:
490 * boxes <-- pointer to set of boxes
491 *
492 * returns:
493 * pointer to global box numbers array
494 *---------------------------------------------------------------------------*/
495
496 const cs_gnum_t *
fvm_box_set_get_g_num(fvm_box_set_t * boxes)497 fvm_box_set_get_g_num(fvm_box_set_t *boxes)
498 {
499 assert(boxes != NULL);
500
501 return boxes->g_num;
502 }
503
504 /*----------------------------------------------------------------------------
505 * Build a Morton_index to get a well-balanced distribution of the boxes.
506 *
507 * parameters:
508 * boxes <-- pointer to associated fvm_box_set_t structure
509 * distrib <-> pointer to a fvm_box_distrib_t structure
510 * n_leaves <-- number of leaves with weight > 0
511 * leaf_codes <-- Morton code for each leaf
512 * weight <-- number of boxes related to each leaf
513 *---------------------------------------------------------------------------*/
514
515 void
fvm_box_set_build_morton_index(const fvm_box_set_t * boxes,fvm_box_distrib_t * distrib,cs_lnum_t n_leaves,fvm_morton_code_t * leaf_codes,cs_lnum_t * weight)516 fvm_box_set_build_morton_index(const fvm_box_set_t *boxes,
517 fvm_box_distrib_t *distrib,
518 cs_lnum_t n_leaves,
519 fvm_morton_code_t *leaf_codes,
520 cs_lnum_t *weight)
521 {
522 #if defined(HAVE_MPI)
523
524 cs_lnum_t *order = NULL;
525
526 assert(distrib != NULL);
527 assert(distrib->morton_index != NULL);
528
529 BFT_MALLOC(order, n_leaves, cs_lnum_t);
530
531 /* Locally order Morton encoding */
532
533 fvm_morton_local_order(n_leaves,
534 leaf_codes,
535 order);
536
537 /* Compute a Morton index on ranks and return the associated fit */
538
539 if (boxes->comm != MPI_COMM_NULL)
540 distrib->fit = fvm_morton_build_rank_index(boxes->dim,
541 distrib->max_level,
542 n_leaves,
543 leaf_codes,
544 weight,
545 order,
546 distrib->morton_index,
547 boxes->comm);
548 /* Free memory */
549
550 BFT_FREE(order);
551
552 #endif
553 }
554
555 /*----------------------------------------------------------------------------
556 * Redistribute boxes over the ranks according to the Morton index to
557 * assume a better balanced distribution of the boxes.
558 *
559 * parameters:
560 * distrib <-- data structure on box distribution
561 * boxes <-> pointer to the structure to redistribute
562 *---------------------------------------------------------------------------*/
563
564 void
fvm_box_set_redistribute(const fvm_box_distrib_t * distrib,fvm_box_set_t * boxes)565 fvm_box_set_redistribute(const fvm_box_distrib_t *distrib,
566 fvm_box_set_t *boxes)
567 {
568 #if defined(HAVE_MPI)
569
570 /* Sanity checks */
571
572 assert(distrib != NULL);
573 assert(boxes != NULL);
574 assert(distrib->n_ranks > 1);
575
576 const cs_lnum_t stride = boxes->dim * 2;
577
578 size_t n_send = distrib->index[distrib->n_ranks];
579
580 int *dest_rank;
581 BFT_MALLOC(dest_rank, n_send, int);
582
583 cs_gnum_t *send_g_num;
584 BFT_MALLOC(send_g_num, n_send, cs_gnum_t);
585 cs_coord_t *send_extents;
586 BFT_MALLOC(send_extents, n_send*stride, cs_coord_t);
587
588 for (int rank_id = 0; rank_id < distrib->n_ranks; rank_id++) {
589 cs_lnum_t s_id = distrib->index[rank_id];
590 cs_lnum_t e_id = distrib->index[rank_id+1];
591 for (cs_lnum_t i = s_id; i < e_id; i++) {
592 cs_lnum_t box_id = distrib->list[i];
593 dest_rank[i] = rank_id;
594 send_g_num[i] = boxes->g_num[box_id];
595 for (cs_lnum_t j = 0; j < stride; j++)
596 send_extents[i*stride + j] = boxes->extents[box_id*stride + j];
597 }
598 }
599
600 BFT_FREE(boxes->g_num);
601 BFT_FREE(boxes->extents);
602
603 cs_all_to_all_t *d = cs_all_to_all_create(n_send,
604 0, /* flags */
605 NULL,
606 dest_rank,
607 boxes->comm);
608
609 /* Exchange global numbers and extents */
610
611 boxes->g_num
612 = cs_all_to_all_copy_array(d,
613 CS_GNUM_TYPE,
614 1,
615 false, /* reverse */
616 send_g_num,
617 NULL);
618
619 boxes->extents
620 = cs_all_to_all_copy_array(d,
621 CS_COORD_TYPE,
622 boxes->dim * 2,
623 false, /* reverse */
624 send_extents,
625 NULL);
626
627 /* Update dimensions */
628
629 boxes->n_boxes = cs_all_to_all_n_elts_dest(d);
630
631 /* Free buffers */
632
633 BFT_FREE(send_extents);
634 BFT_FREE(send_g_num);
635 BFT_FREE(dest_rank);
636
637 cs_all_to_all_destroy(&d);
638
639 #endif /* HAVE_MPI */
640 }
641
642 /*----------------------------------------------------------------------------
643 * Dump a fvm_box_set_t structure.
644 *
645 * parameters:
646 * boxes <-- pointer to the fvm_box_t structure
647 * verbosity <-- verbosity level (0 or 1)
648 *----------------------------------------------------------------------------*/
649
650 void
fvm_box_set_dump(const fvm_box_set_t * boxes,int verbosity)651 fvm_box_set_dump(const fvm_box_set_t *boxes,
652 int verbosity)
653 {
654 cs_lnum_t i;
655
656 const char XYZ[3] = "XYZ";
657
658 if (boxes == NULL)
659 return;
660
661 /* Print basic information */
662
663 if (boxes->dim == 3)
664 bft_printf("\nBox set (3D layout):\n\n"
665 "global min/max on selected faces:\n"
666 " [%7.5e %7.5e %7.5e] --> [%7.5e %7.5e %7.5e]\n",
667 boxes->gmin[0], boxes->gmin[1], boxes->gmin[2],
668 boxes->gmax[0], boxes->gmax[1], boxes->gmax[2]);
669
670 else if (boxes->dim == 2) {
671 bft_printf("\nBox set (2D layout, selected axes [%c, %c]\n\n",
672 XYZ[boxes->dimensions[0]],
673 XYZ[boxes->dimensions[1]]);
674 bft_printf("global min/max on selected faces:\n"
675 " [%7.5e %7.5e] --> [%7.5e %7.5e]\n",
676 boxes->gmin[boxes->dimensions[0]],
677 boxes->gmin[boxes->dimensions[1]],
678 boxes->gmax[boxes->dimensions[0]],
679 boxes->gmax[boxes->dimensions[1]]);
680 }
681
682 else if (boxes->dim == 1) {
683 bft_printf("\nBox set (1D layout, selected axis [%c]\n\n",
684 XYZ[boxes->dimensions[0]]);
685 bft_printf("global min/max on selected faces:\n"
686 " [%7.5e %7.5e] --> [%7.5e %7.5e]\n",
687 boxes->gmin[boxes->dimensions[0]],
688 boxes->gmin[boxes->dimensions[1]],
689 boxes->gmax[boxes->dimensions[0]],
690 boxes->gmax[boxes->dimensions[1]]);
691 }
692 bft_printf_flush();
693
694 /* Print detailed box information */
695
696 if (verbosity < 1)
697 return;
698
699 if (boxes->dim == 3) {
700 for (i = 0; i < boxes->n_boxes; i++) {
701 const cs_coord_t *bmin = boxes->extents + i*6;
702 const cs_coord_t *bmax = boxes->extents + i*6 + 3;
703 bft_printf(" id %8ld, num %9llu: "
704 "[%7.5e %7.5e %7.5e] --> [%7.5e %7.5e %7.5e]\n",
705 (long)i, (unsigned long long)(boxes->g_num[i]),
706 bmin[0], bmin[1], bmin[2],
707 bmax[0], bmax[1], bmax[2]);
708 }
709 }
710
711 else if (boxes->dim == 2) {
712 for (i = 0; i < boxes->n_boxes; i++) {
713 const cs_coord_t *bmin = boxes->extents + i*4;
714 const cs_coord_t *bmax = boxes->extents + i*4 + 2;
715 bft_printf(" id %8ld, num %9llu: "
716 "[%7.5e %7.5e] --> [%7.5e %7.5e]\n",
717 (long)i, (unsigned long long)(boxes->g_num[i]),
718 bmin[0], bmin[1], bmax[0], bmax[1]);
719 }
720 }
721
722 else if (boxes->dim == 1) {
723 for (i = 0; i < boxes->n_boxes; i++) {
724 const cs_coord_t *bmin = boxes->extents + i*2;
725 const cs_coord_t *bmax = boxes->extents + i*2 + 1;
726 bft_printf(" id %8ld, num %9llu: "
727 "[%7.5e] --> [%7.5e]\n",
728 (long)i, (unsigned long long)(boxes->g_num[i]),
729 bmin[0], bmax[0]);
730 }
731 }
732
733 /* Sanity check */
734
735 for (i = 0; i < boxes->n_boxes; i++) {
736 int j;
737 const cs_coord_t *bmin = boxes->extents + boxes->dim*2*i;
738 const cs_coord_t *bmax = boxes->extents + boxes->dim*(2*i + 1);
739 for (j = 0; j < boxes->dim; j++) {
740 if (bmin[j] > bmax[j])
741 bft_error(__FILE__, __LINE__, 0,
742 _("Inconsistent box found (min > max):\n"
743 " global number: %llu\n"
744 " min : %10.4g\n"
745 " max : %10.4g\n"),
746 (unsigned long long)(boxes->g_num[i]), bmin[j], bmax[j]);
747 }
748 }
749
750 }
751
752 #if defined(HAVE_MPI)
753
754 /*----------------------------------------------------------------------------
755 * Create a fvm_box_distrib_t structure.
756 *
757 * parameters:
758 * n_boxes <-- number of boxes
759 * n_g_boxes <-- global number of boxes
760 * max_level <-- max level reached locally in the related tree
761 * comm <-- MPI communicator. on which the distribution takes place
762 *
763 * returns:
764 * a pointer to a new allocated fvm_box_distrib_t structure.
765 *---------------------------------------------------------------------------*/
766
767 fvm_box_distrib_t *
fvm_box_distrib_create(cs_lnum_t n_boxes,cs_gnum_t n_g_boxes,int max_level,MPI_Comm comm)768 fvm_box_distrib_create(cs_lnum_t n_boxes,
769 cs_gnum_t n_g_boxes,
770 int max_level,
771 MPI_Comm comm)
772 {
773 int i, n_ranks, gmax_level;
774
775 fvm_box_distrib_t *new_distrib = NULL;
776
777 if (n_g_boxes == 0)
778 return NULL;
779
780 BFT_MALLOC(new_distrib, 1, fvm_box_distrib_t);
781
782 /* Parallel parameters */
783
784 MPI_Comm_size(comm, &n_ranks);
785
786 new_distrib->n_ranks = n_ranks;
787
788 new_distrib->n_boxes = n_boxes;
789
790 assert(n_ranks > 1);
791
792 BFT_MALLOC(new_distrib->morton_index, n_ranks + 1, fvm_morton_code_t);
793
794 MPI_Allreduce(&max_level, &gmax_level, 1, MPI_INT, MPI_MAX, comm);
795
796 new_distrib->max_level = gmax_level;
797 new_distrib->fit = 999.0;
798
799 BFT_MALLOC(new_distrib->index, n_ranks + 1, cs_lnum_t);
800
801 for (i = 0; i < n_ranks + 1; i++)
802 new_distrib->index[i] = 0;
803
804 new_distrib->list = NULL;
805
806 return new_distrib;
807 }
808
809 /*----------------------------------------------------------------------------
810 * Destroy a fvm_box_distrib_t structure.
811 *
812 * parameters:
813 * distrib <-> pointer to pointer to the structure to destroy
814 *---------------------------------------------------------------------------*/
815
816 void
fvm_box_distrib_destroy(fvm_box_distrib_t ** distrib)817 fvm_box_distrib_destroy(fvm_box_distrib_t **distrib)
818 {
819 if (distrib != NULL) {
820
821 fvm_box_distrib_t *d = *distrib;
822
823 if (d == NULL)
824 return;
825
826 BFT_FREE(d->index);
827 BFT_FREE(d->list);
828 BFT_FREE(d->morton_index);
829
830 BFT_FREE(d);
831 }
832 }
833
834 /*----------------------------------------------------------------------------
835 * Delete redundancies in box distribution
836 *
837 * parameters:
838 * distrib <-> pointer to the fvm_box_distrib_t structure
839 *---------------------------------------------------------------------------*/
840
841 void
fvm_box_distrib_clean(fvm_box_distrib_t * distrib)842 fvm_box_distrib_clean(fvm_box_distrib_t *distrib)
843 {
844 int i, rank;
845
846 cs_lnum_t *counter = NULL, *new_index = NULL;
847
848 BFT_MALLOC(counter, distrib->n_boxes, cs_lnum_t);
849 BFT_MALLOC(new_index, distrib->n_ranks + 1, cs_lnum_t);
850
851 for (i = 0; i < distrib->n_ranks + 1; i++)
852 new_index[i] = 0;
853
854 for (rank = 0; rank < distrib->n_ranks; rank++) {
855
856 cs_lnum_t shift = new_index[rank];
857 cs_lnum_t start = distrib->index[rank];
858 cs_lnum_t end = distrib->index[rank + 1];
859
860 if (end - start > 0) {
861
862 for (i = 0; i < distrib->n_boxes; i++)
863 counter[i] = 0;
864
865 for (i = start; i < end; i++)
866 counter[distrib->list[i]] += 1;
867
868 for (i = 0; i < distrib->n_boxes; i++) {
869
870 if (counter[i] > 0)
871 distrib->list[shift++] = i;
872
873 }
874
875 } /* end if end - start > 0 */
876
877 new_index[rank+1] = shift;
878
879 } /* End of loop on ranks */
880
881 /* Memory management */
882
883 BFT_FREE(distrib->index);
884 BFT_REALLOC(distrib->list, new_index[distrib->n_ranks], cs_lnum_t);
885
886 distrib->index = new_index;
887
888 BFT_FREE(counter);
889 }
890
891 /*----------------------------------------------------------------------------
892 * Display a histogram on leaves associated to the boxes and several
893 * other pieces of information (min, max, ...)
894 *
895 * parameters:
896 * distrib <-- pointer to the fvm_box_distrib_t structure
897 * comm <-- associated MPI communicator
898 *---------------------------------------------------------------------------*/
899
900 void
fvm_box_distrib_dump_statistics(const fvm_box_distrib_t * distrib,MPI_Comm comm)901 fvm_box_distrib_dump_statistics(const fvm_box_distrib_t *distrib,
902 MPI_Comm comm)
903 {
904 cs_lnum_t i;
905
906 int n_ranks = 0;
907 cs_lnum_t n_quantiles = 5;
908 cs_lnum_t quantile_start[6];
909 cs_lnum_t n_boxes[5];
910
911 /* Sanity checks */
912
913 assert(distrib != NULL);
914 assert(distrib->index != NULL);
915
916 _get_distrib_statistics(distrib,
917 &n_quantiles,
918 quantile_start,
919 n_boxes,
920 NULL,
921 &n_ranks,
922 comm);
923
924 bft_printf("\n"
925 "- Box distribution statistics -\n\n");
926
927 bft_printf(" Distribution imbalance: %10.4g\n",
928 distrib->fit);
929 bft_printf(" Number of ranks in distribution: %8d\n\n",
930 n_ranks);
931
932 /* Print histogram to show the distribution of boxes */
933
934 if (n_quantiles > 0) {
935
936 for (i = 0; i < n_quantiles - 1; i++)
937 bft_printf(" %3ld : [ %10ld ; %10ld [ = %10ld\n",
938 (long)i+1, (long)quantile_start[i], (long)quantile_start[i+1],
939 (long)n_boxes[i]);
940
941 i = n_quantiles -1;
942 bft_printf(" %3ld : [ %10ld ; %10ld ] = %10ld\n",
943 (long)i+1, (long)quantile_start[i], (long)quantile_start[i+1] - 1,
944 (long)n_boxes[i]);
945
946 }
947 bft_printf_flush();
948 }
949
950 #endif /* defined(HAVE_MPI) */
951
952 /*---------------------------------------------------------------------------*/
953
954 END_C_DECLS
955