1 /*============================================================================
2 * Functions dealing with ghost cells
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 <stdio.h>
34 #include <string.h>
35 #include <assert.h>
36
37 /*----------------------------------------------------------------------------
38 * Local headers
39 *----------------------------------------------------------------------------*/
40
41 #include "bft_mem.h"
42 #include "bft_error.h"
43 #include "bft_printf.h"
44
45 #include "fvm_periodicity.h"
46
47 #include "cs_base.h"
48 #include "cs_interface.h"
49 #include "cs_mesh.h"
50 #include "cs_order.h"
51 #include "cs_halo.h"
52
53 /*----------------------------------------------------------------------------
54 * Header for the current file
55 *----------------------------------------------------------------------------*/
56
57 #include "cs_mesh_halo.h"
58
59 /*----------------------------------------------------------------------------*/
60
61 BEGIN_C_DECLS
62
63 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
64
65 /*============================================================================
66 * Local structure and type definitions
67 *============================================================================*/
68
69 typedef struct _vtx_lookup_table {
70
71 cs_lnum_t n_vertices; /* Number of local vertices in the problem domain */
72 cs_lnum_t n_transforms; /* Number of transformations */
73 cs_lnum_t n_interfaces; /* Number of interfaces */
74 cs_lnum_t n_categories; /* Number of possible categories
75 = n_interfaces * (n_transforms + 1)
76 (1 category for purely parallel elements) */
77
78 cs_lnum_t *if_ranks; /* List of ranks */
79 cs_lnum_t *rank_ids; /* list of rank ids */
80
81 cs_lnum_t *index; /* index on table (size = n_vertices + 1) */
82
83 cs_lnum_t *rank_list; /* list of ranks on which vertices are linked */
84 cs_lnum_t *type_list; /* list of type (purelly parallel (=0) or number
85 of the periodicity) featuring a vertex. This
86 list is only allocated when n_perio > 0 */
87
88 } vtx_lookup_table_t;
89
90
91 /*=============================================================================
92 * Local Macro definitions
93 *============================================================================*/
94
95 /*============================================================================
96 * Global static variables
97 *============================================================================*/
98
99 /*=============================================================================
100 * Private function definitions
101 *============================================================================*/
102
103 /*---------------------------------------------------------------------------
104 * Fill a look up table structure without periodicity
105 *
106 * The look up list takes the distant rank value with which a local vertex
107 * is linked.
108 *
109 * parameters:
110 * vtx_lookup <-- pointer to a vtx_lookup_table_t structure
111 * ifs <-- pointer to a fvm_interface_set_t structure
112 *---------------------------------------------------------------------------*/
113
114 static void
_fill_vtx_lookup(vtx_lookup_table_t * vtx_lookup,const cs_interface_set_t * ifs)115 _fill_vtx_lookup(vtx_lookup_table_t *vtx_lookup,
116 const cs_interface_set_t *ifs)
117 {
118 cs_lnum_t i, vtx_id, rank_id, shift;
119
120 cs_lnum_t *counter = NULL;
121
122 const cs_lnum_t n_interfaces = cs_interface_set_size(ifs);
123
124 BFT_MALLOC(counter, vtx_lookup->n_vertices, cs_lnum_t);
125
126 for (i = 0; i < vtx_lookup->n_vertices; i++)
127 counter[i] = 0;
128
129 for (rank_id = 0; rank_id < n_interfaces; rank_id++) {
130
131 const cs_interface_t *interface = cs_interface_set_get(ifs, rank_id);
132 const cs_lnum_t interface_size = cs_interface_size(interface);
133 const cs_lnum_t *elt_id = cs_interface_get_elt_ids(interface);
134
135 for (i = 0; i < interface_size; i++) { /* Only parallel vertices */
136
137 vtx_id = elt_id[i];
138 shift = vtx_lookup->index[vtx_id] + counter[vtx_id];
139
140 vtx_lookup->rank_list[shift] = vtx_lookup->rank_ids[rank_id];
141 counter[vtx_id] += 1;
142
143 }
144
145 } /* End of loop on ranks */
146
147 BFT_FREE(counter);
148
149 }
150
151 /*---------------------------------------------------------------------------
152 * Fill a look up table structure with periodicity.
153 *
154 * The rank_list takes the distant rank value with which a local vertex
155 * is linked.
156 * The type_list takes the id number associated to the kind of link between
157 * local and distant element. If the link is purely parallel, list gets value
158 * 0 otherwise the list gets the number of the transformation (+ or -)
159 *
160 * parameters:
161 * vtx_look_up <-- pointer to a vtx_lookup_table_t structure
162 * ifs <-- pointer to a fvm_interface_set_t structure
163 *---------------------------------------------------------------------------*/
164
165 static void
_fill_vtx_lookup_with_perio(vtx_lookup_table_t * vtx_lookup,const cs_interface_set_t * ifs)166 _fill_vtx_lookup_with_perio(vtx_lookup_table_t *vtx_lookup,
167 const cs_interface_set_t *ifs)
168 {
169 cs_lnum_t i, tr_id, vtx_id, rank_id, shift;
170
171 cs_lnum_t *counter = NULL;
172
173 const fvm_periodicity_t *periodicity = cs_interface_set_periodicity(ifs);
174 const cs_lnum_t n_interfaces = cs_interface_set_size(ifs);
175 const cs_lnum_t n_transforms = fvm_periodicity_get_n_transforms(periodicity);
176
177 assert(n_transforms > 0);
178
179 BFT_MALLOC(counter, vtx_lookup->n_vertices, cs_lnum_t);
180
181 for (i = 0; i < vtx_lookup->n_vertices; i++)
182 counter[i] = 0;
183
184 for (rank_id = 0; rank_id < n_interfaces; rank_id++) {
185
186 const cs_interface_t *interface
187 = cs_interface_set_get(ifs, vtx_lookup->rank_ids[rank_id]);
188 const cs_lnum_t tr_index_size = cs_interface_get_tr_index_size(interface);
189 const cs_lnum_t *tr_index = cs_interface_get_tr_index(interface);
190 const cs_lnum_t *vtx_ids = cs_interface_get_elt_ids(interface);
191
192 assert(n_transforms + 2 == tr_index_size);
193 assert(tr_index != NULL);
194
195 for (i = tr_index[0]; i < tr_index[1]; i++) { /* Only parallel vertices */
196
197 vtx_id = vtx_ids[i];
198 shift = vtx_lookup->index[vtx_id] + counter[vtx_id];
199
200 vtx_lookup->rank_list[shift] = rank_id;
201 vtx_lookup->type_list[shift] = 0;
202 counter[vtx_id] += 1;
203
204 }
205
206 for (tr_id = 0; tr_id < n_transforms; tr_id++) {
207
208 for (i = tr_index[tr_id + 1]; i < tr_index[tr_id + 2]; i++) {
209
210 vtx_id = vtx_ids[i];
211 shift = vtx_lookup->index[vtx_id] + counter[vtx_id];
212 vtx_lookup->rank_list[shift] = rank_id;
213 vtx_lookup->type_list[shift] = tr_id + 1;
214 counter[vtx_id] += 1;
215
216 }
217
218 } /* End of loop on transformations */
219
220 } /* End of loop on ranks */
221
222 BFT_FREE(counter);
223
224 }
225
226 /*---------------------------------------------------------------------------
227 * Create a vtx_look_up_table_t structure
228 *
229 * parameters:
230 * n_vertices <-- number of vertices of the table.
231 * ifs <-- pointer to a fvm_interface_set_t structure
232 *
233 * returns:
234 * A pointer to the created vtx_lookup_table_t structure
235 *---------------------------------------------------------------------------*/
236
237 static vtx_lookup_table_t *
_vtx_lookup_create(cs_lnum_t n_vertices,const cs_interface_set_t * ifs)238 _vtx_lookup_create(cs_lnum_t n_vertices,
239 const cs_interface_set_t *ifs)
240 {
241 cs_lnum_t i, rank_id, tmp_id, interface_size;
242
243 cs_lnum_t loc_rank_id = -1;
244 vtx_lookup_table_t *vtx_lookup = NULL;
245
246 const cs_interface_t *interface = NULL;
247 const cs_lnum_t *vtx_ids = NULL;
248 const fvm_periodicity_t *periodicity = cs_interface_set_periodicity(ifs);
249 const cs_lnum_t n_transforms = fvm_periodicity_get_n_transforms(periodicity);
250 const cs_lnum_t n_interfaces = cs_interface_set_size(ifs);
251
252 BFT_MALLOC(vtx_lookup, 1, vtx_lookup_table_t);
253
254 vtx_lookup->n_vertices = n_vertices;
255 vtx_lookup->n_interfaces = n_interfaces;
256 vtx_lookup->n_transforms = n_transforms;
257 vtx_lookup->n_categories = (n_transforms + 1)*n_interfaces;
258
259 BFT_MALLOC(vtx_lookup->index, n_vertices + 1, cs_lnum_t);
260 BFT_MALLOC(vtx_lookup->if_ranks, n_interfaces, cs_lnum_t);
261 BFT_MALLOC(vtx_lookup->rank_ids, n_interfaces, cs_lnum_t);
262
263 for (i = 0; i < n_vertices + 1; i++)
264 vtx_lookup->index[i] = 0;
265
266 /* Check if cs_glob_rank_id belongs to the interface set in order to
267 arrange if_ranks with local rank at first place */
268
269 for (rank_id = 0; rank_id < n_interfaces; rank_id++) {
270
271 interface = cs_interface_set_get(ifs, rank_id);
272 vtx_lookup->if_ranks[rank_id] = cs_interface_rank(interface);
273 vtx_lookup->rank_ids[rank_id] = rank_id;
274
275 if (cs_glob_rank_id == cs_interface_rank(interface))
276 loc_rank_id = rank_id;
277
278 } /* End of loop on if_ranks */
279
280 /* Define vtx_lookup->if_ranks in right order */
281
282 if (loc_rank_id > 0) {
283
284 tmp_id = vtx_lookup->if_ranks[loc_rank_id];
285 vtx_lookup->if_ranks[loc_rank_id] = vtx_lookup->if_ranks[0];
286 vtx_lookup->if_ranks[0] = tmp_id;
287
288 vtx_lookup->rank_ids[0] = loc_rank_id;
289 vtx_lookup->rank_ids[loc_rank_id] = 0;
290
291 /* Order by increasing numbers */
292
293 if (n_interfaces > 2) {
294
295 cs_lnum_t *order = NULL;
296 cs_lnum_t *buffer = NULL;
297 cs_lnum_t *_rank_ids = NULL;
298
299 assert(sizeof(cs_lnum_t) == sizeof(cs_lnum_t));
300
301 BFT_MALLOC(order, n_interfaces - 1, cs_lnum_t);
302 BFT_MALLOC(buffer, n_interfaces - 1, cs_lnum_t);
303 BFT_MALLOC(_rank_ids, n_interfaces , cs_lnum_t);
304
305 _rank_ids[0] = vtx_lookup->rank_ids[0];
306 for (i = 1; i < n_interfaces; i++) {
307 buffer[i-1] = vtx_lookup->if_ranks[i];
308 _rank_ids[i] = vtx_lookup->rank_ids[i];
309 }
310
311 cs_order_lnum_allocated(NULL,
312 buffer,
313 order,
314 n_interfaces-1);
315
316 for (i = 0; i < n_interfaces - 1; i++) {
317 vtx_lookup->if_ranks[i+1] = buffer[order[i]];
318 vtx_lookup->rank_ids[i+1] = _rank_ids[order[i] + 1];
319 }
320
321 BFT_FREE(buffer);
322 BFT_FREE(order);
323 BFT_FREE(_rank_ids);
324
325 } /* End of ordering ranks */
326
327 } /* If rank order has to be changed */
328
329 /* First loop to create index */
330
331 for (rank_id = 0; rank_id < n_interfaces; rank_id++) {
332
333 interface = cs_interface_set_get(ifs, rank_id);
334 interface_size = cs_interface_size(interface);
335 vtx_ids = cs_interface_get_elt_ids(interface);
336
337 for (i = 0; i < interface_size; i++)
338 vtx_lookup->index[vtx_ids[i] + 1] += 1;
339
340 } /* End of loop on if_ranks */
341
342 /* Create index and allocate buffers */
343
344 for (i = 0; i < n_vertices; i++)
345 vtx_lookup->index[i+1] += vtx_lookup->index[i];
346
347 BFT_MALLOC(vtx_lookup->rank_list, vtx_lookup->index[n_vertices], cs_lnum_t);
348
349 /* Second loop to fill table(s) */
350
351 if (n_transforms > 0) {
352
353 BFT_MALLOC(vtx_lookup->type_list, vtx_lookup->index[n_vertices], cs_lnum_t);
354 _fill_vtx_lookup_with_perio(vtx_lookup, ifs);
355
356 }
357 else {
358
359 vtx_lookup->type_list = NULL;
360 _fill_vtx_lookup(vtx_lookup, ifs);
361
362 }
363
364 return vtx_lookup;
365 }
366
367 /*---------------------------------------------------------------------------
368 * Destroy a vtx_lookup structure.
369 *
370 * parameters:
371 * vtx_lookup <-- pointer to a vtx_lookup_table_t structure
372 *---------------------------------------------------------------------------*/
373
374 static void
_vtx_lookup_destroy(vtx_lookup_table_t ** vtx_lookup)375 _vtx_lookup_destroy(vtx_lookup_table_t **vtx_lookup)
376 {
377 vtx_lookup_table_t *vl = *vtx_lookup;
378
379 BFT_FREE(vl->if_ranks);
380 BFT_FREE(vl->rank_ids);
381 BFT_FREE(vl->index);
382 BFT_FREE(vl->rank_list);
383
384 if (vl->type_list != NULL)
385 BFT_FREE(vl->type_list);
386
387 BFT_FREE(*vtx_lookup);
388 }
389
390 /*---------------------------------------------------------------------------
391 * Set checker for this vertex_id according to vtx_lookup features.
392 *
393 * parameters:
394 * vtx_id <-- vertex id to deal with
395 * vtx_checker <-> put a tag in the implied categories
396 * vtx_lookup <-- pointer to a vtx_lookup_table_t structure
397 *---------------------------------------------------------------------------*/
398
399 static void
_update_vtx_checker(cs_lnum_t vtx_id,cs_lnum_t * vtx_checker,vtx_lookup_table_t * vtx_lookup)400 _update_vtx_checker(cs_lnum_t vtx_id,
401 cs_lnum_t *vtx_checker,
402 vtx_lookup_table_t *vtx_lookup)
403 {
404 cs_lnum_t i, rank_id, type;
405
406 const cs_lnum_t n_interfaces = vtx_lookup->n_interfaces;
407 const cs_lnum_t n_transforms = vtx_lookup->n_transforms;
408
409 for (i = vtx_lookup->index[vtx_id];
410 i < vtx_lookup->index[vtx_id + 1]; i++) {
411
412 rank_id = vtx_lookup->rank_list[i];
413
414 if (n_transforms == 0) /* purely parallel */
415
416 vtx_checker[rank_id] += 1;
417
418 else { /* n_perio > 0 */
419
420 type = vtx_lookup->type_list[i];
421 vtx_checker[type*n_interfaces + rank_id] += 1;
422
423 }
424
425 } /* End of loop on vtx_lookup */
426
427 }
428
429 /*---------------------------------------------------------------------------
430 * Update the halo type of each face checker element according to vtx_checker
431 * values for this face.
432 *
433 * parameters:
434 * n_face_vertices <-- number of vertices defining a face
435 * n_categories <-- size of vtx_checker and face_checker
436 * vtx_checker <-- number of vertices for each categories for the face
437 * face_checker --> halo_type value in each categories
438 *
439 * returns:
440 * maximum halo_type value in face_checker
441 *---------------------------------------------------------------------------*/
442
443 static cs_halo_type_t
_update_face_checker(cs_lnum_t n_face_vertices,cs_lnum_t n_categories,cs_lnum_t * vtx_checker,cs_halo_type_t * face_checker)444 _update_face_checker(cs_lnum_t n_face_vertices,
445 cs_lnum_t n_categories,
446 cs_lnum_t *vtx_checker,
447 cs_halo_type_t *face_checker)
448 {
449 cs_lnum_t i;
450
451 cs_halo_type_t ret_type = CS_HALO_N_TYPES;
452
453 for (i = 0; i < n_categories; i++) {
454
455 if (vtx_checker[i] == n_face_vertices) { /* => STANDARD HALO */
456
457 face_checker[i] = CS_HALO_STANDARD;
458 ret_type = CS_HALO_STANDARD;
459
460 }
461 else {
462
463 if (vtx_checker[i] > 0) { /* => EXTENDED HALO */
464
465 if (ret_type == CS_HALO_N_TYPES)
466 ret_type = CS_HALO_EXTENDED;
467
468 if (face_checker[i] == CS_HALO_N_TYPES)
469 face_checker[i] = CS_HALO_EXTENDED;
470
471 }
472
473 }
474
475 } /* End of loop on categories */
476
477 return ret_type;
478 }
479
480 /*---------------------------------------------------------------------------
481 * Update counter on number of elements in each category of halos
482 * according to face_checker values.
483 *
484 * parameters:
485 * mesh <-- pointer to a mesh structure
486 * face_checker <-- halo type in each categories for cell's faces
487 *---------------------------------------------------------------------------*/
488
489 static void
_count_halo_elements(cs_mesh_t * mesh,cs_halo_type_t * face_checker)490 _count_halo_elements(cs_mesh_t *mesh,
491 cs_halo_type_t *face_checker)
492 {
493 cs_lnum_t type_id, rank_id;
494
495 const cs_lnum_t n_transforms = mesh->n_transforms;
496 const cs_halo_t *halo = mesh->halo;
497 const cs_lnum_t n_c_domains = halo->n_c_domains;
498 const cs_lnum_t stride = 4*n_c_domains;
499
500 for (type_id = 0; type_id < n_transforms + 1; type_id++) {
501
502 for (rank_id = 0; rank_id < n_c_domains; rank_id++) {
503
504 if (face_checker[type_id*n_c_domains + rank_id]
505 == CS_HALO_STANDARD) {
506
507 halo->send_index[2*rank_id + 1] += 1;
508
509 if (type_id > 0) /* periodic elements */
510 halo->send_perio_lst[stride*(type_id-1) + 4*rank_id + 1] += 1;
511
512 } /* STANDARD HALO */
513
514 else if (face_checker[type_id*n_c_domains + rank_id]
515 == CS_HALO_EXTENDED) {
516
517 if (mesh->halo_type == CS_HALO_EXTENDED) {
518
519 halo->send_index[2*rank_id + 2] += 1;
520
521 if (type_id > 0) /* periodic elements */
522 halo->send_perio_lst[stride*(type_id-1) + 4*rank_id + 3] += 1;
523
524 }
525
526 } /* EXTENDED HALO */
527
528 } /* End of loop on ranks */
529
530 } /* End of loop on categories */
531
532 }
533
534 /*---------------------------------------------------------------------------
535 * Build halo indexes
536 *
537 * parameters:
538 * mesh <-> pointer to a mesh structure
539 *---------------------------------------------------------------------------*/
540
541 static void
_build_halo_index(cs_mesh_t * mesh)542 _build_halo_index(cs_mesh_t *mesh)
543 {
544 cs_lnum_t i, rank_id;
545 cs_lnum_t buffer_size, n_halo_elts, n_per_halo_elts;
546
547 cs_halo_t *halo = mesh->halo;
548
549 const cs_lnum_t n_c_domains = halo->n_c_domains;
550 const cs_lnum_t n_init_perio = mesh->n_init_perio;
551 const cs_lnum_t stride = 4*n_c_domains;
552 const cs_lnum_t n_transforms = mesh->n_transforms;
553 const cs_lnum_t local_rank = (cs_glob_rank_id == -1) ? 0:cs_glob_rank_id;
554
555 buffer_size = 0;
556 for (i = 0; i < 2*n_c_domains; i++)
557 buffer_size += halo->send_index[i+1];
558
559 cs_alloc_mode_t halo_alloc_mode = cs_halo_get_buffer_alloc_mode();
560
561 CS_MALLOC_HD(halo->send_list, buffer_size, cs_lnum_t, halo_alloc_mode);
562
563 /* Define parallel and periodic index */
564
565 for (rank_id = 0; rank_id < n_c_domains; rank_id++) {
566
567 /* Standard halo */
568 /* ------------- */
569
570 /* Count number of periodic elements in standard halo for this rank */
571
572 n_per_halo_elts = 0;
573 for (i = 0; i < n_transforms; i++)
574 n_per_halo_elts += halo->send_perio_lst[stride*i + 4*rank_id + 1];
575
576 /* Define index */
577
578 n_halo_elts = halo->send_index[2*rank_id+1];
579 halo->send_index[2*rank_id+1] += halo->send_index[2*rank_id];
580
581 assert(n_halo_elts >= n_per_halo_elts);
582
583 /* Fill perio_lst buffer */
584
585 if (n_init_perio > 0) {
586
587 if (halo->c_domain_rank[rank_id] == local_rank)
588 halo->send_perio_lst[4*rank_id] = halo->send_index[2*rank_id];
589
590 else
591 halo->send_perio_lst[4*rank_id] =
592 halo->send_index[2*rank_id] + (n_halo_elts - n_per_halo_elts);
593
594 for (i = 0; i < n_transforms - 1; i++)
595 halo->send_perio_lst[stride*(i+1) + 4*rank_id]
596 = halo->send_perio_lst[stride*i + 4*rank_id]
597 + halo->send_perio_lst[stride*i + 4*rank_id + 1];
598
599 } /* Test if n_perio > 0 */
600
601 /* Extended halo */
602 /* ------------- */
603
604 n_per_halo_elts = 0;
605 for (i = 0; i < n_transforms; i++)
606 n_per_halo_elts += halo->send_perio_lst[stride*i + 4*rank_id+3];
607
608 n_halo_elts = halo->send_index[2*rank_id+2];
609 halo->send_index[2*rank_id+2] += halo->send_index[2*rank_id+1];
610
611 assert(n_halo_elts >= n_per_halo_elts);
612
613 if (n_init_perio > 0) {
614
615 if (halo->c_domain_rank[rank_id] == local_rank)
616 halo->send_perio_lst[4*rank_id+2] = halo->send_index[2*rank_id+1];
617
618 else
619 halo->send_perio_lst[4*rank_id+2] =
620 halo->send_index[2*rank_id+1] + (n_halo_elts - n_per_halo_elts);
621
622 for (i = 0; i < n_transforms - 1; i++)
623 halo->send_perio_lst[stride*(i+1) + 4*rank_id + 2]
624 = halo->send_perio_lst[stride*i + 4*rank_id + 2]
625 + halo->send_perio_lst[stride*i + 4*rank_id + 3];
626
627 } /* Test if n_perio > 0 */
628
629 } /* End of loop on c_domain_rank */
630
631 }
632
633 /*---------------------------------------------------------------------------
634 * Fill ghost cells list (member of cs_halo_t structure)
635 *
636 * parameters:
637 * mesh <-- pointer to a mesh structure.
638 * face_checker <-- halo type of each face of the cell.
639 * cell_id <-- numbering of the treated cell.
640 * counter <-> counter on each categories.
641 *---------------------------------------------------------------------------*/
642
643 static void
_add_halo_elements(cs_mesh_t * mesh,cs_halo_type_t * face_checker,cs_lnum_t cell_id,cs_lnum_t * counter)644 _add_halo_elements(cs_mesh_t *mesh,
645 cs_halo_type_t *face_checker,
646 cs_lnum_t cell_id,
647 cs_lnum_t *counter)
648 {
649 cs_lnum_t i, type_id, shift, c_shift;
650
651 cs_halo_t *halo = mesh->halo;
652
653 const cs_lnum_t n_transforms = mesh->n_transforms;
654 const cs_lnum_t n_c_domains = halo->n_c_domains;
655
656 /* How is organized counter:
657
658 -------------------------------------------------
659 Paral: | | | | | | | | | | | | |
660 -------------------------------------------------
661 std ext std ext std ext std ext std ext std ext
662 _______ _______ _______ _______ _______ _______
663 rank0 rank1 rank2 rank3 rank4 rank5
664
665 -------------------------------------------------
666 P1: | | | | | | | | | | | | |
667 -------------------------------------------------
668 std ext std ext std ext std ext std ext std ext
669 _______ _______ _______ _______ _______ _______
670 rank0 rank1 rank2 rank3 rank4 rank5
671
672 -------------------------------------------------
673 P-1: | | | | | | | | | | | | |
674 -------------------------------------------------
675 std ext std ext std ext std ext std ext std ext
676 _______ _______ _______ _______ _______ _______
677 rank0 rank1 rank2 rank3 rank4 rank5
678
679 etc...
680
681 */
682
683 for (type_id = 0; type_id < n_transforms + 1; type_id++) {
684
685 for (i = 0; i < n_c_domains; i++) {
686
687 if (face_checker[type_id*n_c_domains + i] == CS_HALO_STANDARD) {
688
689 c_shift = 2*n_c_domains*type_id + 2*i;
690
691 if (type_id == 0)
692 shift = halo->send_index[2*i] + counter[c_shift];
693 else
694 shift = halo->send_perio_lst[4*n_c_domains*(type_id-1) + 4*i]
695 + counter[c_shift];
696
697 halo->send_list[shift] = cell_id;
698 counter[c_shift] += 1;
699
700 }
701 else if (face_checker[type_id*n_c_domains + i] == CS_HALO_EXTENDED) {
702
703 if (mesh->halo_type == CS_HALO_EXTENDED) {
704
705 c_shift = 2*n_c_domains*type_id + 2*i + 1;
706
707 if (type_id == 0)
708 shift = halo->send_index[2*i+1] + counter[c_shift];
709 else
710 shift = halo->send_perio_lst[4*n_c_domains*(type_id-1) + 4*i + 2]
711 + counter[c_shift];
712
713 halo->send_list[shift] = cell_id;
714 counter[c_shift] += 1;
715
716 } /* End of extended halo treatment */
717
718 }
719
720 } /* End of loop on interfaces */
721
722 } /* End of loop on categories */
723
724 }
725
726 /*---------------------------------------------------------------------------
727 * Test if loop has to be continued according halo type and face number.
728 *
729 * parameters:
730 * mesh <-- pointer to mesh structure
731 * face_id <-- interior face id
732 *---------------------------------------------------------------------------*/
733
734 static bool
_test_loop_continues(cs_mesh_t * mesh,cs_lnum_t face_id)735 _test_loop_continues(cs_mesh_t *mesh,
736 cs_lnum_t face_id)
737 {
738 bool choice = false;
739
740 /* Face has to be an internal face */
741
742 if (mesh->halo_type == CS_HALO_STANDARD) {
743
744 if ( mesh->i_face_cells[face_id][0] < 0
745 || mesh->i_face_cells[face_id][1] < 0)
746 choice = true;
747 else
748 choice = false;
749
750 }
751 else {
752
753 assert(mesh->halo_type == CS_HALO_EXTENDED);
754 choice = true;
755
756 }
757
758 return choice;
759 }
760
761 /*---------------------------------------------------------------------------
762 * Define the elements of send_halo structure.
763 *
764 * Two main loops. First one for counting number of elements and create index.
765 * Second one for filling the ghost cells list.
766 *
767 * parameters:
768 * mesh <-- pointer to cs_mesh_t structure
769 * vertex_ifs <-- pointer to fvm_vertex_ifs_t structure
770 * gcell_faces_idx <-- "ghost cell -> faces" connectivity index
771 * gcell_faces_lst <-- "ghost cell -> faces" connectivity list
772 *---------------------------------------------------------------------------*/
773
774 static void
_fill_send_halo(cs_mesh_t * mesh,const cs_interface_set_t * vertex_ifs,cs_lnum_t * gcell_faces_idx,cs_lnum_t * gcell_faces_lst)775 _fill_send_halo(cs_mesh_t *mesh,
776 const cs_interface_set_t *vertex_ifs,
777 cs_lnum_t *gcell_faces_idx,
778 cs_lnum_t *gcell_faces_lst)
779 {
780 cs_lnum_t i, cell_id, i_fac, i_vtx;
781 cs_lnum_t fac_id, vtx_id;
782 cs_lnum_t n_face_vertices;
783
784 cs_halo_type_t type_tag = CS_HALO_N_TYPES;
785 cs_halo_type_t face_type = CS_HALO_N_TYPES;
786 cs_halo_type_t cell_type = CS_HALO_N_TYPES;
787 cs_lnum_t n_categories = 0;
788 cs_halo_t *halo = mesh->halo;
789 vtx_lookup_table_t *vtx_lookup = NULL;
790 cs_halo_type_t *cell_tag = NULL;
791 cs_halo_type_t *face_checker = NULL;
792 cs_lnum_t *vtx_checker = NULL;
793 cs_lnum_t *counter = NULL;
794
795 const cs_lnum_t n_cells = mesh->n_cells;
796 const cs_lnum_t n_vertices = mesh->n_vertices;
797 const cs_lnum_t *fac_vtx_idx = mesh->i_face_vtx_idx;
798 const cs_lnum_t *fac_vtx_lst = mesh->i_face_vtx_lst;
799
800 /* We should have the faces -> vertices connectivity to continue */
801
802 if (mesh->i_face_vtx_lst == NULL)
803 return;
804
805 /* Create a lookup table to accelerate search in
806 cs_interface_set structure */
807
808 vtx_lookup = _vtx_lookup_create(n_vertices, vertex_ifs);
809
810 n_categories = vtx_lookup->n_categories;
811
812 /* How is organized vtx_checker and face_checker:
813
814 -------------------------
815 Paral: | | | | |
816 -------------------------
817 rank0 rank1 rank2 rank3
818
819 -------------------------
820 P1: | | | | |
821 -------------------------
822 rank0 rank1 rank2 rank3
823
824 -------------------------
825 P-1: | | | | |
826 -------------------------
827 rank0 rank1 rank2 rank3
828
829 */
830
831 if (gcell_faces_idx != NULL && gcell_faces_lst != NULL) {
832
833 BFT_MALLOC(vtx_checker, n_categories, cs_lnum_t);
834 BFT_MALLOC(face_checker, n_categories, cs_halo_type_t);
835 BFT_MALLOC(cell_tag, n_cells, cs_halo_type_t);
836
837 /* First loop to create index and allocate cell_list */
838
839 for (cell_id = 0; cell_id < n_cells; cell_id++) {
840
841 /* Default initialization */
842
843 cell_type = CS_HALO_N_TYPES;
844
845 for (i = 0; i < n_categories; i++)
846 face_checker[i] = CS_HALO_N_TYPES;
847
848 /* Loop on faces of the cell */
849
850 for (i_fac = gcell_faces_idx[cell_id];
851 i_fac < gcell_faces_idx[cell_id + 1];
852 i_fac++) {
853
854 fac_id = gcell_faces_lst[i_fac];
855
856 if (_test_loop_continues(mesh, fac_id) == true) {
857
858 n_face_vertices = fac_vtx_idx[fac_id + 1] - fac_vtx_idx[fac_id];
859
860 /* Initialize checker */
861
862 for (i = 0; i < n_categories; i++)
863 vtx_checker[i] = 0;
864
865 /* Loop on vertices of the face */
866
867 for (i_vtx = fac_vtx_idx[fac_id];
868 i_vtx < fac_vtx_idx[fac_id + 1];
869 i_vtx++) {
870
871 vtx_id = fac_vtx_lst[i_vtx];
872
873 _update_vtx_checker(vtx_id, vtx_checker, vtx_lookup);
874
875 } /* End of loop on vertices */
876
877 face_type = _update_face_checker(n_face_vertices,
878 n_categories,
879 vtx_checker,
880 face_checker);
881
882 cell_type = CS_MIN(face_type, cell_type);
883
884 } /* If the face is an internal face */
885
886 } /* End of loop on faces */
887
888 _count_halo_elements(mesh, face_checker);
889
890 cell_tag[cell_id] = cell_type;
891
892 } /* End of loop on cells */
893
894 /* Build halo index */
895
896 _build_halo_index(mesh);
897
898 /* Initialize counter */
899
900 BFT_MALLOC(counter, 2*n_categories, cs_lnum_t);
901
902 for (i = 0; i < 2*n_categories; i++)
903 counter[i] = 0;
904
905 /* Second loop to build halo->ghost_cells */
906
907 for (cell_id = 0; cell_id < n_cells; cell_id++) {
908
909 if (mesh->halo_type == CS_HALO_STANDARD)
910 type_tag = CS_HALO_EXTENDED;
911 else if (mesh->halo_type == CS_HALO_EXTENDED)
912 type_tag = CS_HALO_N_TYPES;
913
914 if (cell_tag[cell_id] < type_tag) {
915
916 for (i = 0; i < n_categories; i++)
917 face_checker[i] = CS_HALO_N_TYPES;
918
919 /* Loop on faces of the cell */
920
921 for (i_fac = gcell_faces_idx[cell_id];
922 i_fac < gcell_faces_idx[cell_id+1];
923 i_fac++) {
924
925 /* Initialize checker */
926
927 for (i = 0; i < n_categories; i++)
928 vtx_checker[i] = 0;
929
930 fac_id = gcell_faces_lst[i_fac];
931
932 if (_test_loop_continues(mesh, fac_id) == true) {
933
934 /* Loop on vertices of the face */
935
936 n_face_vertices = fac_vtx_idx[fac_id + 1] - fac_vtx_idx[fac_id];
937
938 for (i_vtx = fac_vtx_idx[fac_id];
939 i_vtx < fac_vtx_idx[fac_id + 1];
940 i_vtx++) {
941
942 vtx_id = fac_vtx_lst[i_vtx];
943
944 _update_vtx_checker(vtx_id,
945 vtx_checker,
946 vtx_lookup);
947
948 } /* End of loop on vertices */
949
950 face_type = _update_face_checker(n_face_vertices,
951 n_categories,
952 vtx_checker,
953 face_checker);
954
955 } /* If face is an internal face */
956
957 } /* End of loop on faces */
958
959 _add_halo_elements(mesh,
960 face_checker,
961 cell_id,
962 counter);
963
964 } /* End of test on cell_list */
965
966 } /* End of loop on cells */
967
968 /* Free memory */
969
970 BFT_FREE(vtx_checker);
971 BFT_FREE(face_checker);
972 BFT_FREE(counter);
973 BFT_FREE(cell_tag);
974
975 } /* End if gcell_face_idx and gcell_faces_lst != NULL */
976
977 /* Destroy the lookup table structure */
978
979 _vtx_lookup_destroy(&vtx_lookup);
980
981 /* Complete halo definition */
982
983 halo->n_send_elts[CS_HALO_STANDARD] = 0;
984 halo->n_send_elts[CS_HALO_EXTENDED] = 0;
985
986 for (i = 0; i < halo->n_c_domains; i++) {
987
988 halo->n_send_elts[CS_HALO_STANDARD] += halo->send_index[2*i+1]
989 - halo->send_index[2*i];
990 halo->n_send_elts[CS_HALO_EXTENDED] += halo->send_index[2*i+2]
991 - halo->send_index[2*i+1];
992
993 }
994
995 halo->n_send_elts[CS_HALO_EXTENDED] += halo->n_send_elts[CS_HALO_STANDARD];
996 }
997
998 /*---------------------------------------------------------------------------
999 * Define a buffer on vertices where vertex belonging to the vertex_ifs
1000 * are tagged with 1 else 0.
1001 *
1002 * parameters:
1003 * n_vertices <-- size of the buffer
1004 * vertex_ifs <-- pointer to a cs_interface_set_t structure
1005 * p_vertex_tag <-> pointer to the tagged buffer
1006 *---------------------------------------------------------------------------*/
1007
1008 static void
_get_vertex_tag(cs_lnum_t n_vertices,const cs_interface_set_t * vertex_ifs,cs_lnum_t * p_vertex_tag[])1009 _get_vertex_tag(cs_lnum_t n_vertices,
1010 const cs_interface_set_t *vertex_ifs,
1011 cs_lnum_t *p_vertex_tag[])
1012 {
1013 cs_lnum_t i, j, rank_id;
1014
1015 cs_lnum_t *vertex_tag = NULL;
1016
1017 const int ifs_size = cs_interface_set_size(vertex_ifs);
1018
1019 BFT_MALLOC(vertex_tag, n_vertices, cs_lnum_t);
1020
1021 for (i = 0; i < n_vertices; i++)
1022 vertex_tag[i] = 0;
1023
1024 for (rank_id = 0; rank_id < ifs_size; rank_id++) {
1025
1026 const cs_interface_t *interface = cs_interface_set_get(vertex_ifs, rank_id);
1027 const cs_lnum_t *vtx_ids = cs_interface_get_elt_ids(interface);
1028 const cs_lnum_t if_size = cs_interface_size(interface);
1029
1030 for (j = 0; j < if_size; j++)
1031 vertex_tag[vtx_ids[j]] = 1;
1032
1033 } /* End of loop on ranks */
1034
1035 *p_vertex_tag = vertex_tag;
1036
1037 }
1038
1039 /*---------------------------------------------------------------------------
1040 * Compute the number of purely parallel ghost cells for a specific rank.
1041 *
1042 * parameters:
1043 * mesh <-- pointer to cs_mesh_t structure
1044 * rank <-- rank on which we want to know the number of purely
1045 * parallel elements.
1046 * type <-- standard or extended
1047 * index <-- index on halo's elements
1048 * perio_lst <-- periodic details on halo
1049 *
1050 * returns:
1051 * Number of purely parallel elements in the halo.
1052 *---------------------------------------------------------------------------*/
1053
1054 static cs_lnum_t
_get_n_par_ghost_cells(const cs_mesh_t * mesh,cs_lnum_t rank,cs_halo_type_t type,const cs_lnum_t index[],const cs_lnum_t perio_lst[])1055 _get_n_par_ghost_cells(const cs_mesh_t *mesh,
1056 cs_lnum_t rank,
1057 cs_halo_type_t type,
1058 const cs_lnum_t index[],
1059 const cs_lnum_t perio_lst[])
1060 {
1061 cs_lnum_t i;
1062 cs_lnum_t n_per_gcells = 0, n_par_gcells = 0;
1063
1064 const cs_lnum_t n_transforms = mesh->n_transforms;
1065 const cs_lnum_t n_c_domains = mesh->halo->n_c_domains;
1066
1067 if (type == CS_HALO_STANDARD) {
1068
1069 for (i = 0; i < n_transforms; i++)
1070 n_per_gcells += perio_lst[4*rank+1 + 4*n_c_domains*i];
1071
1072 n_par_gcells = index[2*rank+1] - index[2*rank];
1073 n_par_gcells -= n_per_gcells;
1074
1075 }
1076 else if (type == CS_HALO_EXTENDED) {
1077
1078 for (i = 0; i < n_transforms; i++)
1079 n_per_gcells += perio_lst[4*rank+3 + 4*n_c_domains*i];
1080
1081 n_par_gcells = index[2*rank+2] - index[2*rank+1];
1082 n_par_gcells -= n_per_gcells;
1083
1084 }
1085
1086 return n_par_gcells;
1087 }
1088
1089 /*---------------------------------------------------------------------------
1090 * Exchange number and list of cells constituting send_halo structure for each
1091 * frontier ranks. Fill the halo structure from these data.
1092 *
1093 * parameters:
1094 * mesh <-- pointer to a cs_mesh_t structure
1095 *---------------------------------------------------------------------------*/
1096
1097 static void
_fill_halo(cs_mesh_t * mesh)1098 _fill_halo(cs_mesh_t *mesh)
1099 {
1100 cs_lnum_t rank_id, i, j;
1101 cs_lnum_t shift;
1102
1103 #if defined(HAVE_MPI)
1104 MPI_Request _request[128];
1105 MPI_Request *request = _request;
1106 MPI_Status _status[128];
1107 MPI_Status *status = _status;
1108 #endif
1109
1110 int request_count = 0;
1111
1112 cs_lnum_t *count = NULL;
1113
1114 cs_halo_t *halo = mesh->halo;
1115
1116 const cs_lnum_t n_c_domains = halo->n_c_domains;
1117 const cs_lnum_t n_transforms = mesh->n_transforms;
1118 const cs_lnum_t local_rank = (cs_glob_rank_id == -1) ? 0:cs_glob_rank_id;
1119
1120 #if defined(HAVE_MPI)
1121 if (halo->n_c_domains*2 > 128) {
1122 BFT_MALLOC(request, halo->n_c_domains*2, MPI_Request);
1123 BFT_MALLOC(status, halo->n_c_domains*2, MPI_Status);
1124 }
1125 #endif
1126
1127 /* Build index */
1128 /* ----------- */
1129 /* Receive data from distant ranks */
1130
1131 for (rank_id = 0; rank_id < n_c_domains; rank_id++) {
1132
1133 if (halo->c_domain_rank[rank_id] != local_rank) {
1134
1135 #if defined(HAVE_MPI)
1136 MPI_Irecv(&(halo->index[2*rank_id+1]), 2, CS_MPI_LNUM,
1137 halo->c_domain_rank[rank_id],
1138 halo->c_domain_rank[rank_id],
1139 cs_glob_mpi_comm,
1140 &(request[request_count++]));
1141 #endif
1142
1143 }
1144
1145 } /* End of loop on ranks */
1146
1147 /* We wait for receiving all messages */
1148
1149 #if defined(HAVE_MPI)
1150 if (mesh->n_domains > 1)
1151 MPI_Barrier(cs_glob_mpi_comm);
1152 #endif
1153
1154 BFT_MALLOC(count, 2*n_c_domains, cs_lnum_t);
1155
1156 /* Send data to distant ranks */
1157
1158 for (rank_id = 0; rank_id < n_c_domains; rank_id++) {
1159
1160 shift = 2*rank_id;
1161 count[shift] = halo->send_index[2*rank_id+1]
1162 - halo->send_index[2*rank_id];
1163 count[shift+1] = halo->send_index[2*rank_id+2]
1164 - halo->send_index[2*rank_id+1];
1165
1166 if (halo->c_domain_rank[rank_id] != local_rank) {
1167
1168 #if defined(HAVE_MPI)
1169 MPI_Isend(&(count[shift]), 2, CS_MPI_LNUM,
1170 halo->c_domain_rank[rank_id],
1171 local_rank,
1172 cs_glob_mpi_comm,
1173 &(request[request_count++]));
1174 #endif
1175
1176 }
1177 else {
1178
1179 halo->index[shift+1] = count[shift];
1180 halo->index[shift+2] = count[shift+1];
1181
1182 }
1183
1184 } /* End of loop on ranks */
1185
1186 /* Wait for all exchanges being done */
1187
1188 #if defined(HAVE_MPI)
1189 if (mesh->n_domains > 1)
1190 MPI_Waitall(request_count, request, status);
1191 #endif
1192 request_count = 0;
1193
1194 BFT_FREE(count);
1195
1196 /* Build index */
1197 /*-------------*/
1198
1199 for (i = 0; i < 2*n_c_domains; i++)
1200 halo->index[i+1] += halo->index[i];
1201
1202 /* Exchange number of elements for each periodicity and for each rank.
1203 Then build halo->perio_lst */
1204
1205 if (mesh->n_init_perio > 0) {
1206
1207 cs_lnum_t n_elts;
1208 cs_lnum_t *exchange_buffer = NULL;
1209
1210 /* n_transforms periodicities to deal with and for each sub-periodicity
1211 2 data. One for standard halo and the other one for extended halo */
1212
1213 const cs_lnum_t n_elts_to_exchange = 2 * n_transforms;
1214
1215 BFT_MALLOC(exchange_buffer, 4*n_transforms, cs_lnum_t);
1216
1217 for (rank_id = 0; rank_id < n_c_domains; rank_id++) {
1218
1219 if (halo->c_domain_rank[rank_id] != local_rank) {
1220
1221 /* Fill buffer to send */
1222
1223 for (i = 0; i < n_transforms; i++) {
1224 shift = 4*n_c_domains*i + 4*rank_id;
1225 for (j = 0; j < 2; j++)
1226 exchange_buffer[2*i+j] = halo->send_perio_lst[shift + 2*j + 1];
1227 }
1228
1229 #if defined(HAVE_MPI)
1230 MPI_Sendrecv(&(exchange_buffer[0]), n_elts_to_exchange, CS_MPI_LNUM,
1231 halo->c_domain_rank[rank_id], local_rank,
1232 &(exchange_buffer[n_elts_to_exchange]), n_elts_to_exchange,
1233 CS_MPI_LNUM,
1234 halo->c_domain_rank[rank_id], halo->c_domain_rank[rank_id],
1235 cs_glob_mpi_comm, status);
1236 #endif
1237
1238 /* Put received elements in the periodic structure */
1239
1240 for (i = 0; i < n_transforms; i++) {
1241 shift = 4*n_c_domains*i + 4*rank_id;
1242 for (j = 0; j < 2; j++)
1243 halo->perio_lst[shift + 2*j + 1] =
1244 exchange_buffer[n_elts_to_exchange + 2*i + j];
1245 }
1246
1247 } /* rank != local_rank */
1248
1249 else {
1250
1251 for (i = 0; i < n_transforms; i++) {
1252
1253 shift = 4*n_c_domains*i + 4*rank_id;
1254 for (j = 0; j < 2; j++)
1255 halo->perio_lst[shift + 2*j + 1] =
1256 halo->send_perio_lst[shift + 2*j + 1];
1257
1258 } /* End of loop on periodicities */
1259
1260 } /* local_rank == rank */
1261
1262 } /* End of loop on communicating ranks */
1263
1264 BFT_FREE(exchange_buffer);
1265
1266 /* Build index values for perio_lst */
1267
1268 for (rank_id = 0; rank_id < n_c_domains; rank_id++) {
1269
1270 /* Build index for standard ghost cells */
1271
1272 n_elts = _get_n_par_ghost_cells(mesh,
1273 rank_id,
1274 CS_HALO_STANDARD,
1275 halo->index,
1276 halo->perio_lst);
1277
1278 halo->perio_lst[4*rank_id] = halo->index[2*rank_id] + n_elts;
1279
1280 for (i = 0; i < n_transforms - 1; i++) {
1281 shift = 4*n_c_domains*i + 4*rank_id;
1282 halo->perio_lst[4*n_c_domains + shift] =
1283 halo->perio_lst[shift] + halo->perio_lst[shift + 1];
1284 }
1285
1286 /* Build index for extended ghost cells */
1287
1288 n_elts = _get_n_par_ghost_cells(mesh,
1289 rank_id,
1290 CS_HALO_EXTENDED,
1291 halo->index,
1292 halo->perio_lst);
1293
1294 halo->perio_lst[4*rank_id+2] = halo->index[2*rank_id+1] + n_elts;
1295
1296 for (i = 0; i < n_transforms - 1; i++) {
1297 shift = 4*n_c_domains*i + 4*rank_id + 2;
1298 halo->perio_lst[4*n_c_domains + shift] =
1299 halo->perio_lst[shift] + halo->perio_lst[shift + 1];
1300 }
1301
1302 } /* End of loop on communicating ranks */
1303
1304 } /* End if n_perio > 0 */
1305
1306 #if defined(HAVE_MPI)
1307 if (request != _request) {
1308 BFT_FREE(request);
1309 BFT_FREE(status);
1310 }
1311 #endif
1312
1313 halo->n_elts[CS_HALO_STANDARD] = 0;
1314 halo->n_elts[CS_HALO_EXTENDED] = 0;
1315
1316 for (i = 0; i < n_c_domains; i++) {
1317
1318 halo->n_elts[CS_HALO_STANDARD] += halo->index[2*i+1] - halo->index[2*i];
1319 halo->n_elts[CS_HALO_EXTENDED] += halo->index[2*i+2] - halo->index[2*i+1];
1320
1321 }
1322
1323 halo->n_elts[CS_HALO_EXTENDED] += halo->n_elts[CS_HALO_STANDARD];
1324 }
1325
1326 /*---------------------------------------------------------------------------
1327 * Compute maximum list buffer size.
1328 *
1329 * This is done to avoid a reallocation for each rank and transformation.
1330 *
1331 * parameters:
1332 * ifs <-- pointer to a cs_interface_set_t structure
1333 *
1334 * returns:
1335 * max buffer size
1336 *---------------------------------------------------------------------------*/
1337
1338 static cs_lnum_t
_get_list_buffer_size(const cs_interface_set_t * ifs)1339 _get_list_buffer_size(const cs_interface_set_t *ifs)
1340 {
1341 cs_lnum_t i, j, tr_index_size;
1342
1343 cs_lnum_t max_lst_size = 0;
1344
1345 const cs_interface_t *interface = NULL;
1346 const cs_lnum_t *tr_index = NULL;
1347 const cs_lnum_t ifs_size = cs_interface_set_size(ifs);
1348
1349 if (ifs == NULL)
1350 return max_lst_size;
1351
1352 for (i = 0; i < ifs_size; i++) {
1353
1354 interface = cs_interface_set_get(ifs, i);
1355 tr_index = cs_interface_get_tr_index(interface);
1356 tr_index_size = cs_interface_get_tr_index_size(interface) - 1;
1357
1358 if (tr_index != NULL)
1359 for (j = 0; j < tr_index_size; j++)
1360 max_lst_size = CS_MAX(max_lst_size, tr_index[j+1] - tr_index[j]);
1361 else
1362 max_lst_size = CS_MAX(max_lst_size,
1363 (cs_lnum_t)cs_interface_size(interface));
1364
1365 } /* End of loop on interfaces */
1366
1367 return max_lst_size;
1368 }
1369
1370 /*---------------------------------------------------------------------------
1371 * Define an index on vertices belonging to this interface for this rank
1372 * and this transformation.
1373 *
1374 * parameters:
1375 * ifs <-- pointer to cs_interface_set_t structure
1376 * rank_id <-- rank number to work with
1377 * tr_id <-- transformation id to work with
1378 * vtx_interface_idx <-> index on vertices which match the criterions
1379 *---------------------------------------------------------------------------*/
1380
1381 static void
_define_vtx_interface_idx(const cs_interface_set_t * ifs,cs_lnum_t rank_id,cs_lnum_t tr_id,cs_lnum_t n_vertices,cs_lnum_t * vtx_interface_idx)1382 _define_vtx_interface_idx(const cs_interface_set_t *ifs,
1383 cs_lnum_t rank_id,
1384 cs_lnum_t tr_id,
1385 cs_lnum_t n_vertices,
1386 cs_lnum_t *vtx_interface_idx)
1387 {
1388 cs_lnum_t i, j, id;
1389
1390 /* Initialize index */
1391
1392 for (i = 0; i < n_vertices + 1; i++)
1393 vtx_interface_idx[i] = 0;
1394
1395 for (id = 0; id < cs_interface_set_size(ifs); id++) {
1396
1397 const cs_interface_t *interface = cs_interface_set_get(ifs, id);
1398
1399 if (rank_id == cs_interface_rank(interface)) {
1400
1401 const cs_lnum_t *tr_index = cs_interface_get_tr_index(interface);
1402 const cs_lnum_t *vtx_ids = cs_interface_get_elt_ids(interface);
1403
1404 if (tr_index == NULL) { /* purely parallel elements */
1405
1406 for (j = 0; j < (cs_lnum_t)cs_interface_size(interface); j++)
1407 vtx_interface_idx[vtx_ids[j] + 1] += 1;
1408
1409 }
1410 else {
1411
1412 for (j = tr_index[tr_id]; j < tr_index[tr_id+1]; j++)
1413 vtx_interface_idx[vtx_ids[j] + 1] += 1;
1414
1415 }
1416
1417 /* Create index */
1418
1419 for (j = 0; j < n_vertices; j++)
1420 vtx_interface_idx[j+1] += vtx_interface_idx[j];
1421
1422 break;
1423
1424 }
1425
1426 } /* End of loop on interfaces */
1427
1428 }
1429
1430 /*---------------------------------------------------------------------------
1431 * Define an index on vertices belonging the interface for this rank
1432 * and this perio.
1433 * Fill the dist_id_lst which is the list of distant ids associated
1434 * to local vertices (same rank and same transformation).
1435 *
1436 * parameters:
1437 * ifs <-- pointer to cs_interface_set_t structure
1438 * rank_id <-- rank number to work with
1439 * tr_id <-- transformation id to work with
1440 * n_vertices <-- number of vertices
1441 * dist_id_lst <-> list of distant vertex numbers matching criteria
1442 * vtx_interface_idx <-> index on vertices matching criteria
1443 *---------------------------------------------------------------------------*/
1444
1445 static void
_define_dist_id_lst(const cs_interface_set_t * ifs,int rank_id,int tr_id,cs_lnum_t n_vertices,cs_lnum_t * dist_id_lst,cs_lnum_t * vtx_interface_idx)1446 _define_dist_id_lst(const cs_interface_set_t *ifs,
1447 int rank_id,
1448 int tr_id,
1449 cs_lnum_t n_vertices,
1450 cs_lnum_t *dist_id_lst,
1451 cs_lnum_t *vtx_interface_idx)
1452 {
1453 cs_lnum_t i, j, id, shift;
1454
1455 /* Initialize index */
1456
1457 for (i = 0; i < n_vertices + 1; i++)
1458 vtx_interface_idx[i] = 0;
1459
1460 for (id = 0; id < cs_interface_set_size(ifs); id++) {
1461
1462 const cs_interface_t *interface = cs_interface_set_get(ifs, id);
1463
1464 if (rank_id == cs_interface_rank(interface)) {
1465
1466 const cs_lnum_t *tr_index = cs_interface_get_tr_index(interface);
1467 const cs_lnum_t *vtx_ids = cs_interface_get_elt_ids(interface);
1468 const cs_lnum_t *dist_vtx_ids = cs_interface_get_match_ids(interface);
1469
1470 if (tr_index == NULL)
1471 for (j = 0; j < (cs_lnum_t)cs_interface_size(interface); j++)
1472 vtx_interface_idx[vtx_ids[j] + 1] += 1;
1473
1474 else
1475 for (j = tr_index[tr_id]; j < tr_index[tr_id+1]; j++)
1476 vtx_interface_idx[vtx_ids[j] + 1] += 1;
1477
1478 /* Create index */
1479
1480 for (j = 0; j < n_vertices; j++)
1481 vtx_interface_idx[j+1] += vtx_interface_idx[j];
1482
1483 /* There must be only one distant vertex id per vtx_id when
1484 we handle a specific rank and a specific periodicity.
1485 So, we don't need a counter to fill dist_id_lst */
1486
1487 if (tr_index == NULL) {
1488
1489 for (j = 0; j < (cs_lnum_t)cs_interface_size(interface); j++) {
1490 shift = vtx_interface_idx[vtx_ids[j]];
1491 dist_id_lst[shift] = dist_vtx_ids[j];
1492 }
1493
1494
1495 }
1496 else {
1497
1498 for (j = tr_index[tr_id]; j < tr_index[tr_id+1]; j++) {
1499 shift = vtx_interface_idx[vtx_ids[j]];
1500 dist_id_lst[shift] = dist_vtx_ids[j];
1501 }
1502
1503 }
1504
1505 break;
1506
1507 }
1508
1509 } /* End of loop on interfaces */
1510
1511 }
1512
1513 /*---------------------------------------------------------------------------
1514 * Compute the start and end index in ghost cells list for the send_halo
1515 * elements according to its rank, its periodicity and its type.
1516 *
1517 * parameters:
1518 * mesh <-- pointer to cs_mesh_t structure
1519 * index <-- index on halo's elements
1520 * perio_lst <-- periodic details on halo
1521 * rank_id <-- rank number to work with
1522 * tr_id <-- transformation id to work with
1523 * type_id <-- standard or extended
1524 * p_start_idx --> pointer on start index
1525 * p_end_idx --> pointer on end index
1526 *---------------------------------------------------------------------------*/
1527
1528 static void
_get_start_end_idx(const cs_mesh_t * mesh,const cs_lnum_t * index,const cs_lnum_t * perio_lst,int rank_id,cs_lnum_t tr_id,cs_lnum_t type_id,cs_lnum_t * p_start_idx,cs_lnum_t * p_end_idx)1529 _get_start_end_idx(const cs_mesh_t *mesh,
1530 const cs_lnum_t *index,
1531 const cs_lnum_t *perio_lst,
1532 int rank_id,
1533 cs_lnum_t tr_id,
1534 cs_lnum_t type_id,
1535 cs_lnum_t *p_start_idx,
1536 cs_lnum_t *p_end_idx)
1537 {
1538 cs_lnum_t i, n_par_gcells, n_per_gcells;
1539 cs_lnum_t start_idx = -1, end_idx = -1;
1540
1541 const cs_lnum_t n_c_domains = mesh->halo->n_c_domains;
1542
1543 if (tr_id == 0) { /* Purelly parallel elements */
1544
1545 if (type_id == 0) { /* STANDARD HALO */
1546
1547 n_par_gcells = _get_n_par_ghost_cells(mesh,
1548 rank_id,
1549 CS_HALO_STANDARD,
1550 index,
1551 perio_lst);
1552
1553 start_idx = index[2*rank_id];
1554 end_idx = start_idx + n_par_gcells;
1555
1556 }
1557
1558 if (type_id == 1) { /* EXTENDED HALO */
1559
1560 n_par_gcells = _get_n_par_ghost_cells(mesh,
1561 rank_id,
1562 CS_HALO_EXTENDED,
1563 index,
1564 perio_lst);
1565
1566 start_idx = index[2*rank_id+1];
1567 end_idx = start_idx + n_par_gcells;
1568
1569 }
1570
1571 }
1572 else { /* Periodic elements */
1573
1574 i = tr_id - 1;
1575
1576 if (type_id == 0) { /* STANDARD HALO */
1577
1578 n_per_gcells = perio_lst[4*rank_id + 1 + 4*n_c_domains*i];
1579 start_idx = perio_lst[4*rank_id + 4*n_c_domains*i];
1580 end_idx = start_idx + n_per_gcells;
1581
1582 }
1583
1584 if (type_id == 1) { /* EXTENDED HALO */
1585
1586 n_per_gcells = perio_lst[4*rank_id + 3 + 4*n_c_domains*i];
1587 start_idx = perio_lst[4*rank_id + 2 + 4*n_c_domains*i];
1588 end_idx = start_idx + n_per_gcells;
1589
1590 }
1591
1592 } /* Parallel or periodic elements */
1593
1594 *p_start_idx = start_idx;
1595 *p_end_idx = end_idx;
1596
1597 }
1598
1599 /*---------------------------------------------------------------------------
1600 * Compute the size of the "ghost cell to distant vertices" connectivity for
1601 * send_halo elements.
1602 * This will be use to define "ghost cell to distant vertices" index.
1603 *
1604 * parameters:
1605 * mesh <-- pointer to cs_mesh_t structure
1606 * ifs <-- pointer to cs_interface_set_t structure
1607 * rank_id <-- rank number to work with
1608 * tr_id <-- transformation id to work with
1609 * cell_faces_idx <-- "cell -> faces" connectivity index
1610 * cell_faces_lst <-- "cell -> faces" connectivity list
1611 * vtx_interface_idx <-> index on vertices which match the criterions
1612 * vtx_tag <-> tag array on vertices
1613 * gcell_dist_vtx_idx <-> "ghost cell -> distant vertices" connectivity index
1614 *---------------------------------------------------------------------------*/
1615
1616 static void
_count_send_gcell_to_dist_vtx_connect(cs_mesh_t * mesh,const cs_interface_set_t * ifs,int rank_id,int tr_id,cs_lnum_t * cell_faces_idx,cs_lnum_t * cell_faces_lst,cs_lnum_t * vtx_interface_idx,cs_lnum_t * vtx_tag,cs_lnum_t * gcell_dist_vtx_idx)1617 _count_send_gcell_to_dist_vtx_connect(cs_mesh_t *mesh,
1618 const cs_interface_set_t *ifs,
1619 int rank_id,
1620 int tr_id,
1621 cs_lnum_t *cell_faces_idx,
1622 cs_lnum_t *cell_faces_lst,
1623 cs_lnum_t *vtx_interface_idx,
1624 cs_lnum_t *vtx_tag,
1625 cs_lnum_t *gcell_dist_vtx_idx)
1626 {
1627 cs_lnum_t id, cell_id, i_fac, i_vtx, i_loop;
1628 cs_lnum_t start_idx, end_idx, vtx_id, fac_id;
1629
1630 cs_lnum_t n_loops = 0;
1631 cs_lnum_t n_added_vertices = 0;
1632
1633 cs_halo_t *halo = mesh->halo;
1634
1635 const cs_lnum_t n_vertices = mesh->n_vertices;
1636 const cs_lnum_t *fac_vtx_idx = mesh->i_face_vtx_idx;
1637 const cs_lnum_t *fac_vtx_lst = mesh->i_face_vtx_lst;
1638
1639 const char err_corresp[]
1640 = N_("Repeated inconsistency in the halo construction.\n"
1641 "Several local points have the same distant correspondant;\n"
1642 "this is probably a side effect of a poor quality joining\n"
1643 "or periodicity, or of a partitioning with periodicity issue.\n"
1644 "Coordinates of the first impacted point: [%12.5e, %12.5e %12.5e].");
1645
1646 _define_vtx_interface_idx(ifs,
1647 halo->c_domain_rank[rank_id],
1648 tr_id,
1649 n_vertices,
1650 vtx_interface_idx);
1651
1652 /* Count size of the connectivity and define index */
1653
1654 if (mesh->halo_type == CS_HALO_STANDARD)
1655 n_loops = 1;
1656 else if (mesh->halo_type == CS_HALO_EXTENDED)
1657 n_loops = 2;
1658
1659 for (i_loop = 0; i_loop < n_loops; i_loop++) {
1660
1661 /* Define start and end idx */
1662
1663 _get_start_end_idx(mesh,
1664 halo->send_index,
1665 halo->send_perio_lst,
1666 rank_id,
1667 tr_id,
1668 i_loop,
1669 &start_idx,
1670 &end_idx);
1671
1672 for (id = start_idx; id < end_idx; id++) {
1673
1674 cell_id = halo->send_list[id];
1675
1676 for (i_fac = cell_faces_idx[cell_id];
1677 i_fac < cell_faces_idx[cell_id+1];
1678 i_fac++) {
1679
1680 fac_id = cell_faces_lst[i_fac];
1681
1682 if (_test_loop_continues(mesh, fac_id) == true) {
1683
1684 /* Loop on vertices of the face */
1685
1686 for (i_vtx = fac_vtx_idx[fac_id];
1687 i_vtx < fac_vtx_idx[fac_id+1];
1688 i_vtx++) {
1689
1690 vtx_id = fac_vtx_lst[i_vtx];
1691
1692 /* If vertex is on the interface for this rank and this
1693 transformation */
1694
1695 n_added_vertices = vtx_interface_idx[vtx_id+1]
1696 - vtx_interface_idx[vtx_id];
1697
1698 if (n_added_vertices > 0) {
1699
1700 if (n_added_vertices > 1) {
1701
1702 if (cs_glob_n_ranks > 1)
1703 bft_printf("fac_num: %ld (%llu)\n"
1704 "vtx_num: %ld (%llu) - n_added: %ld\n",
1705 (long)fac_id+1,
1706 (unsigned long long)(mesh->global_i_face_num[fac_id]),
1707 (long)vtx_id+1,
1708 (unsigned long long)(mesh->global_vtx_num[vtx_id]),
1709 (long)n_added_vertices);
1710 else
1711 bft_printf("fac_num: %ld\n"
1712 "vtx_num: %ld - n_added: %ld\n",
1713 (long)fac_id+1, (long)vtx_id+1,
1714 (long)n_added_vertices);
1715 bft_printf_flush();
1716
1717 bft_error(__FILE__, __LINE__, 0, _(err_corresp),
1718 mesh->vtx_coord[vtx_id*3],
1719 mesh->vtx_coord[vtx_id*3+1],
1720 mesh->vtx_coord[vtx_id*3+2]);
1721
1722 }
1723
1724 /* Add this vertex if not already checked */
1725
1726 if (vtx_tag[vtx_id] != id) {
1727 vtx_tag[vtx_id] = id;
1728 gcell_dist_vtx_idx[id+1] += n_added_vertices;
1729 }
1730
1731 }
1732
1733 } /* End of loop on vertices */
1734
1735 } /* Treatment only for implied faces */
1736
1737 } /* End of loop on faces */
1738
1739 } /* End of loop on ghost cells */
1740
1741 } /* Loop on STANDARD and EXTENDED halo */
1742
1743 }
1744
1745 /*---------------------------------------------------------------------------
1746 * Fill the "ghost cells to distant vertices" connectivity.
1747 *
1748 * parameters:
1749 * mesh <-- pointer to cs_mesh_t structure
1750 * ifs <-- pointer to cs_interface_set_t structure
1751 * rank_id <-- rank number to work with
1752 * tr_id <-- transformation id to work with
1753 * cell_faces_idx <-- "cell -> faces" connectivity index
1754 * cell_faces_lst <-- "cell -> faces" connectivity list
1755 * vtx_interface_idx <-> index on vertices matching criteria
1756 * dist_id_lst <-> list of distant vertex numbers matching criteria
1757 * counter <-> counter on vertices
1758 * vtx_tag <-> tag array on vertices
1759 * gcell_dist_vtx_idx <-> "ghost cell -> distant vertices" connectivity index
1760 * gcell_dist_vtx_lst <-> "ghost cell -> distant vertices" connectivity list
1761 *---------------------------------------------------------------------------*/
1762
1763 static void
_fill_send_gcell_to_dist_vtx_connect(cs_mesh_t * mesh,const cs_interface_set_t * ifs,cs_lnum_t rank_id,cs_lnum_t tr_id,cs_lnum_t * cell_faces_idx,cs_lnum_t * cell_faces_lst,cs_lnum_t * vtx_interface_idx,cs_lnum_t * dist_id_lst,cs_lnum_t * counter,cs_lnum_t * vtx_tag,cs_lnum_t * gcell_dist_vtx_idx,cs_lnum_t * gcell_dist_vtx_lst)1764 _fill_send_gcell_to_dist_vtx_connect(cs_mesh_t *mesh,
1765 const cs_interface_set_t *ifs,
1766 cs_lnum_t rank_id,
1767 cs_lnum_t tr_id,
1768 cs_lnum_t *cell_faces_idx,
1769 cs_lnum_t *cell_faces_lst,
1770 cs_lnum_t *vtx_interface_idx,
1771 cs_lnum_t *dist_id_lst,
1772 cs_lnum_t *counter,
1773 cs_lnum_t *vtx_tag,
1774 cs_lnum_t *gcell_dist_vtx_idx,
1775 cs_lnum_t *gcell_dist_vtx_lst)
1776 {
1777 cs_lnum_t i, id, cell_id, i_fac, i_vtx, i_loop;
1778 cs_lnum_t shift, vtx_id, fac_id, start_idx, end_idx;
1779
1780 cs_lnum_t n_loops = 0;
1781 cs_halo_t *halo = mesh->halo;
1782
1783 const cs_lnum_t n_vertices = mesh->n_vertices;
1784 const cs_lnum_t *fac_vtx_idx = mesh->i_face_vtx_idx;
1785 const cs_lnum_t *fac_vtx_lst = mesh->i_face_vtx_lst;
1786
1787 _define_dist_id_lst(ifs,
1788 halo->c_domain_rank[rank_id],
1789 tr_id,
1790 n_vertices,
1791 dist_id_lst,
1792 vtx_interface_idx);
1793
1794 /* Fill the "ghost cells to distant vertices" connectivity */
1795
1796 if (mesh->halo_type == CS_HALO_STANDARD)
1797 n_loops = 1;
1798 else if (mesh->halo_type == CS_HALO_EXTENDED)
1799 n_loops = 2;
1800
1801 for (i_loop = 0; i_loop < n_loops; i_loop++) {
1802
1803 /* Define start and end idx */
1804
1805 _get_start_end_idx(mesh,
1806 halo->send_index,
1807 halo->send_perio_lst,
1808 rank_id,
1809 tr_id,
1810 i_loop,
1811 &start_idx,
1812 &end_idx);
1813
1814 for (id = start_idx; id < end_idx; id++) {
1815
1816 cell_id = halo->send_list[id];
1817
1818 for (i_fac = cell_faces_idx[cell_id];
1819 i_fac < cell_faces_idx[cell_id+1];
1820 i_fac++) {
1821
1822 fac_id = cell_faces_lst[i_fac];
1823
1824 if (_test_loop_continues(mesh, fac_id) == true) {
1825
1826 /* Loop on vertices of the face */
1827
1828 for (i_vtx = fac_vtx_idx[fac_id];
1829 i_vtx < fac_vtx_idx[fac_id + 1];
1830 i_vtx++) {
1831
1832 vtx_id = fac_vtx_lst[i_vtx];
1833
1834 /* If vertex is on the interface for this rank and periodicity */
1835
1836 if (vtx_interface_idx[vtx_id+1] - vtx_interface_idx[vtx_id] > 0) {
1837
1838 /* Add this vertex if nont already checked */
1839
1840 if (vtx_tag[vtx_id] != id) { /* Add this vertex */
1841
1842 vtx_tag[vtx_id] = id;
1843
1844 for (i = vtx_interface_idx[vtx_id];
1845 i < vtx_interface_idx[vtx_id+1];
1846 i++) {
1847
1848 shift = gcell_dist_vtx_idx[id] + counter[id];
1849 gcell_dist_vtx_lst[shift] = dist_id_lst[i];
1850 counter[id] += 1;
1851
1852 }
1853
1854 }
1855
1856 } /* If there is something to fill */
1857
1858 } /* End of loop on vertices */
1859
1860 } /* Treatment only for implied faces */
1861
1862 } /* End of loop on faces */
1863
1864 } /* End of loop on ghost cells */
1865
1866 } /* Loop on STANDARD or EXTENDED halo */
1867
1868 }
1869
1870 /*---------------------------------------------------------------------------
1871 * Create a local "ghost cells -> distant vertices" connectivity for
1872 * send_halo cells.
1873 *
1874 * parameters:
1875 * mesh <-- pointer to cs_mesh_t structure
1876 * interface_set <-- pointer to cs_interface_set_t structure
1877 * cell_faces_idx <-- "cell -> faces" connectivity index
1878 * cell_faces_lst <-- "cell -> faces" connectivity list
1879 * p_gcell_dist_vtx_idx --> "ghost cell -> distant vertices" connect. index
1880 * p_gcell_dist_vtx_lst --> "ghost cell -> distant vertices" connect. list
1881 *---------------------------------------------------------------------------*/
1882
1883 static void
_create_send_gcell_vtx_connect(cs_mesh_t * mesh,cs_interface_set_t * interface_set,cs_lnum_t * cell_faces_idx,cs_lnum_t * cell_faces_lst,cs_lnum_t * p_gcell_dist_vtx_idx[],cs_lnum_t * p_gcell_dist_vtx_lst[])1884 _create_send_gcell_vtx_connect(cs_mesh_t *mesh,
1885 cs_interface_set_t *interface_set,
1886 cs_lnum_t *cell_faces_idx,
1887 cs_lnum_t *cell_faces_lst,
1888 cs_lnum_t *p_gcell_dist_vtx_idx[],
1889 cs_lnum_t *p_gcell_dist_vtx_lst[])
1890 {
1891 cs_lnum_t i, id, rank_id;
1892
1893 cs_lnum_t *gcell_dist_vtx_idx = NULL, *gcell_dist_vtx_lst = NULL;
1894 cs_lnum_t *vtx_interface_idx = NULL;
1895 cs_lnum_t *dist_id_lst = NULL;
1896 cs_lnum_t *vtx_tag = NULL;
1897 cs_lnum_t *counter = NULL;
1898
1899 cs_halo_t *halo = mesh->halo;
1900
1901 const cs_lnum_t max_lst_size = _get_list_buffer_size(interface_set);
1902 const cs_lnum_t n_ghost_cells = halo->n_send_elts[CS_HALO_EXTENDED];
1903 const cs_lnum_t n_vertices = mesh->n_vertices;
1904 const cs_lnum_t n_c_domains = halo->n_c_domains;
1905 const cs_lnum_t tr_index_size = mesh->n_transforms + 1;
1906
1907 /* Allocate and initialize buffers */
1908
1909 BFT_MALLOC(gcell_dist_vtx_idx, n_ghost_cells + 1, cs_lnum_t);
1910 BFT_MALLOC(counter, n_ghost_cells, cs_lnum_t);
1911
1912 gcell_dist_vtx_idx[0] = 0;
1913 for (i = 0; i < n_ghost_cells; i++) {
1914 gcell_dist_vtx_idx[i+1] = 0;
1915 counter[i] = 0;
1916 }
1917
1918 BFT_MALLOC(vtx_tag, n_vertices, cs_lnum_t);
1919
1920 for (i = 0; i < n_vertices; i++)
1921 vtx_tag[i] = -1;
1922
1923 BFT_MALLOC(vtx_interface_idx, n_vertices + 1, cs_lnum_t);
1924 BFT_MALLOC(dist_id_lst, max_lst_size, cs_lnum_t);
1925
1926 for (i = 0; i < max_lst_size; i++)
1927 dist_id_lst[i] = -1;
1928
1929 /* Loop on each rank belonging to send_halo.
1930 Create a vertex to ghost cells connectivity for each rank */
1931
1932 for (rank_id = 0; rank_id < n_c_domains; rank_id++) {
1933
1934 /* Define gcell_dist_vtx_idx */
1935
1936 for (id = 0; id < tr_index_size; id++)
1937 _count_send_gcell_to_dist_vtx_connect(mesh,
1938 interface_set,
1939 rank_id,
1940 id,
1941 cell_faces_idx,
1942 cell_faces_lst,
1943 vtx_interface_idx,
1944 vtx_tag,
1945 gcell_dist_vtx_idx);
1946
1947 } /* End of loop on ranks */
1948
1949 /* Create gcell_dist_vtx_idx */
1950
1951 for (i = 0; i < n_ghost_cells; i++)
1952 gcell_dist_vtx_idx[i+1] += gcell_dist_vtx_idx[i];
1953
1954 BFT_MALLOC(gcell_dist_vtx_lst, gcell_dist_vtx_idx[n_ghost_cells], cs_lnum_t);
1955
1956 for (i = 0; i < n_vertices; i++)
1957 vtx_tag[i] = -1;
1958
1959 cs_interface_set_add_match_ids(interface_set);
1960
1961 for (rank_id = 0; rank_id < n_c_domains; rank_id++) {
1962
1963 for (id = 0; id < tr_index_size; id++)
1964 _fill_send_gcell_to_dist_vtx_connect(mesh,
1965 interface_set,
1966 rank_id,
1967 id,
1968 cell_faces_idx,
1969 cell_faces_lst,
1970 vtx_interface_idx,
1971 dist_id_lst,
1972 counter,
1973 vtx_tag,
1974 gcell_dist_vtx_idx,
1975 gcell_dist_vtx_lst);
1976
1977 } /* End of loop on ranks */
1978
1979 cs_interface_set_free_match_ids(interface_set);
1980
1981 BFT_FREE(counter);
1982 BFT_FREE(vtx_tag);
1983 BFT_FREE(vtx_interface_idx);
1984 BFT_FREE(dist_id_lst);
1985
1986 *p_gcell_dist_vtx_idx = gcell_dist_vtx_idx;
1987 *p_gcell_dist_vtx_lst = gcell_dist_vtx_lst;
1988
1989 }
1990
1991 /*---------------------------------------------------------------------------
1992 * Send "ghost cells to distant_num vertices" connectivity on communicating
1993 * ranks and receive the same kind of connectivity from distant ranks.
1994 *
1995 * parameters:
1996 * mesh <-- pointer to cs_mesh_t structure
1997 * send_gcell_dist_vtx_idx --> "ghost cell -> distant vertices" index
1998 * send_gcell_dist_vtx_lst --> "ghost cell -> distant vertices" list
1999 * p_gcell_dist_vtx_idx <-- "ghost cell -> distant vertices" index
2000 * p_gcell_dist_vtx_lst <-- "ghost cell -> distant vertices" list
2001 *---------------------------------------------------------------------------*/
2002
2003 static void
_exchange_gcell_vtx_connect(cs_mesh_t * mesh,cs_lnum_t * send_gcell_dist_vtx_idx,cs_lnum_t * send_gcell_dist_vtx_lst,cs_lnum_t * p_gcell_dist_vtx_idx[],cs_lnum_t * p_gcell_dist_vtx_lst[])2004 _exchange_gcell_vtx_connect(cs_mesh_t *mesh,
2005 cs_lnum_t *send_gcell_dist_vtx_idx,
2006 cs_lnum_t *send_gcell_dist_vtx_lst,
2007 cs_lnum_t *p_gcell_dist_vtx_idx[],
2008 cs_lnum_t *p_gcell_dist_vtx_lst[])
2009 {
2010 cs_lnum_t i, j, rank_id;
2011 cs_lnum_t send_start_idx, send_end_idx, start_idx, end_idx;
2012 cs_lnum_t n_send_elts, n_recv_elts;
2013
2014 cs_lnum_t send_buffer_size = 0;
2015
2016 cs_lnum_t *send_idx_buffer = NULL;
2017 cs_lnum_t *gcell_dist_vtx_idx = NULL, *gcell_dist_vtx_lst = NULL;
2018 cs_lnum_t *send_buffer = NULL, *recv_buffer = NULL;
2019
2020 cs_halo_t *halo = mesh->halo;
2021
2022 const cs_lnum_t local_rank = (cs_glob_rank_id == -1) ? 0:cs_glob_rank_id;
2023 const cs_lnum_t n_c_domains = halo->n_c_domains;
2024 const cs_lnum_t n_ghost_cells = halo->n_elts[CS_HALO_EXTENDED];
2025
2026 #if defined(HAVE_MPI)
2027 MPI_Status status;
2028 #endif
2029
2030 /* Allocate buffers */
2031
2032 for (rank_id = 0; rank_id < n_c_domains; rank_id++) {
2033 if (halo->c_domain_rank[rank_id] != local_rank) {
2034 n_send_elts = halo->send_index[2*rank_id+2]- halo->send_index[2*rank_id];
2035 send_buffer_size = CS_MAX(send_buffer_size, n_send_elts);
2036 }
2037 }
2038
2039 BFT_MALLOC(send_idx_buffer, send_buffer_size, cs_lnum_t);
2040 BFT_MALLOC(gcell_dist_vtx_idx, n_ghost_cells + 1, cs_lnum_t);
2041
2042 for (i = 0; i < n_ghost_cells + 1; i++)
2043 gcell_dist_vtx_idx[i] = 0;
2044
2045 /* Exchange sizes to define gcell_dist_vtx_idx */
2046
2047 for (rank_id = 0; rank_id < n_c_domains; rank_id++) {
2048
2049 recv_buffer = &(gcell_dist_vtx_idx[1 + halo->index[2*rank_id]]);
2050
2051 if (halo->c_domain_rank[rank_id] != local_rank) {
2052
2053 /* Fill send buffer */
2054
2055 for (i = halo->send_index[2*rank_id], j = 0;
2056 i < halo->send_index[2*rank_id+2]; i++, j++)
2057 send_idx_buffer[j] = ( send_gcell_dist_vtx_idx[i+1]
2058 - send_gcell_dist_vtx_idx[i]);
2059
2060 n_send_elts = halo->send_index[2*rank_id+2] - halo->send_index[2*rank_id];
2061 n_recv_elts = halo->index[2*rank_id+2] - halo->index[2*rank_id];
2062
2063 #if defined(HAVE_MPI)
2064 MPI_Sendrecv(&(send_idx_buffer[0]), n_send_elts, CS_MPI_LNUM,
2065 halo->c_domain_rank[rank_id], local_rank,
2066 &(recv_buffer[0]), n_recv_elts, CS_MPI_LNUM,
2067 halo->c_domain_rank[rank_id], halo->c_domain_rank[rank_id],
2068 cs_glob_mpi_comm, &status);
2069 #endif
2070
2071 } /* If rank != local_rank */
2072
2073 else {
2074
2075 for (i = halo->send_index[2*rank_id], j = 0;
2076 i < halo->send_index[2*rank_id+2]; i++, j++)
2077 recv_buffer[j] = ( send_gcell_dist_vtx_idx[i+1]
2078 - send_gcell_dist_vtx_idx[i]);
2079
2080 } /* rank == local_rank */
2081
2082 } /* End of loop on if_ranks */
2083
2084 BFT_FREE(send_idx_buffer);
2085
2086 /* Define index */
2087
2088 for (i = 0; i < n_ghost_cells; i++)
2089 gcell_dist_vtx_idx[i+1] += gcell_dist_vtx_idx[i];
2090
2091 BFT_MALLOC(gcell_dist_vtx_lst, gcell_dist_vtx_idx[n_ghost_cells], cs_lnum_t);
2092
2093 /* Exchange lists to define gcell_dist_vtx_lst */
2094
2095 for (rank_id = 0; rank_id < n_c_domains; rank_id++) {
2096
2097 /* Exchange conectivity list */
2098
2099 send_start_idx = send_gcell_dist_vtx_idx[halo->send_index[2*rank_id]];
2100 send_end_idx = send_gcell_dist_vtx_idx[halo->send_index[2*rank_id+2]];
2101 n_send_elts = send_end_idx - send_start_idx;
2102 send_buffer = &(send_gcell_dist_vtx_lst[send_start_idx]);
2103
2104 start_idx = gcell_dist_vtx_idx[halo->index[2*rank_id]];
2105 end_idx = gcell_dist_vtx_idx[halo->index[2*rank_id+2]];
2106 n_recv_elts = end_idx - start_idx;
2107 recv_buffer = &(gcell_dist_vtx_lst[start_idx]);
2108
2109 if (halo->c_domain_rank[rank_id] != local_rank) {
2110
2111 #if defined(HAVE_MPI)
2112 MPI_Sendrecv(send_buffer, n_send_elts, CS_MPI_LNUM,
2113 halo->c_domain_rank[rank_id], local_rank,
2114 recv_buffer, n_recv_elts, CS_MPI_LNUM,
2115 halo->c_domain_rank[rank_id], halo->c_domain_rank[rank_id],
2116 cs_glob_mpi_comm, &status);
2117 #endif
2118
2119 }
2120 else {
2121
2122 assert(n_recv_elts == n_send_elts);
2123
2124 for (i = 0; i < n_send_elts; i++)
2125 recv_buffer[i] = send_buffer[i];
2126
2127 }
2128
2129 } /* End of loop on ranks */
2130
2131 *p_gcell_dist_vtx_idx = gcell_dist_vtx_idx;
2132 *p_gcell_dist_vtx_lst = gcell_dist_vtx_lst;
2133
2134 }
2135
2136 /*---------------------------------------------------------------------------
2137 * Check mesh->i_face_cells array for ghost cells in standard halo.
2138 *
2139 * parameters:
2140 * mesh <-> pointer to cs_mesh_t structure
2141 *---------------------------------------------------------------------------*/
2142
2143 static void
_check_i_face_cells(cs_mesh_t * mesh)2144 _check_i_face_cells(cs_mesh_t *mesh)
2145 {
2146 int i;
2147
2148 for (i = 0; i < mesh->n_i_faces; i++) {
2149
2150 if (mesh->i_face_cells[i][0] < 0) {
2151 if (mesh->global_i_face_num != NULL)
2152 bft_error(__FILE__, __LINE__, 0,
2153 " Error detected in interior face -> cells connectivity.\n"
2154 " Face %ld (%llu) has an incomplete connectivity.\n"
2155 " Cell1: %ld - Cell2: %ld (%llu)",
2156 (long)i+1, (unsigned long long)(mesh->global_i_face_num[i]),
2157 (long)mesh->i_face_cells[i][0],
2158 (long)mesh->i_face_cells[i][1],
2159 (unsigned long long)(mesh->global_cell_num[mesh->i_face_cells[i][1]]));
2160 else /* Serial run */
2161 bft_error(__FILE__, __LINE__, 0,
2162 " Error detected in interior face -> cells connectivity.\n"
2163 " Face %ld has an incomplete connectivity.\n"
2164 " Cell1: %ld - Cell2: %ld",
2165 (long)i+1, (long)mesh->i_face_cells[i][0],
2166 (long)mesh->i_face_cells[i][1]);
2167 }
2168
2169 if (mesh->i_face_cells[i][1] < 0) {
2170 if (mesh->global_i_face_num != NULL)
2171 bft_error(__FILE__, __LINE__, 0,
2172 " Error detected in interior face -> cells connectivity.\n"
2173 " Face %ld (%llu) has an incomplete connectivity.\n"
2174 " Cell1: %ld (%llu) - Cell2: %ld",
2175 (long)i+1, (unsigned long long)(mesh->global_i_face_num[i]),
2176 (long)mesh->i_face_cells[i][0],
2177 (unsigned long long)(mesh->global_cell_num[mesh->i_face_cells[i][0]]),
2178 (long)mesh->i_face_cells[i][1]);
2179 else
2180 bft_error(__FILE__, __LINE__, 0,
2181 " Error detected in interior face -> cells connectivity.\n"
2182 " Face %ld has an incomplete connectivity.\n"
2183 " Cell1: %ld - Cell2: %ld",
2184 (long)i+1, (long)mesh->i_face_cells[i][0],
2185 (long)mesh->i_face_cells[i][1]);
2186 }
2187
2188 }
2189
2190 }
2191
2192 /*---------------------------------------------------------------------------
2193 * Tag local face by its distant face id and reset previous tag
2194 *
2195 * parameters:
2196 * prev <-- previous start index
2197 * start <-- start index in lnum and dnum
2198 * end <-- end index in lnum and dnum
2199 * l_ids <-- elt_ids array in cs_interface_t struct.
2200 * d_ids <-- match_ids array in cs_interface_t struct.
2201 * l2d_fids <-> tag on local faces
2202 *---------------------------------------------------------------------------*/
2203
2204 static void
_tag_interface_faces(int prev,int start,int end,const cs_lnum_t l_ids[],const cs_lnum_t d_ids[],cs_lnum_t l2d_fids[])2205 _tag_interface_faces(int prev,
2206 int start,
2207 int end,
2208 const cs_lnum_t l_ids[],
2209 const cs_lnum_t d_ids[],
2210 cs_lnum_t l2d_fids[])
2211 {
2212 cs_lnum_t i;
2213
2214 assert(l2d_fids != NULL);
2215
2216 if (prev > -1) /* Reset previous face num. */
2217 for (i = prev; i < start; i++)
2218 l2d_fids[l_ids[i]] = -1;
2219
2220 for (i = start; i < end; i++)
2221 l2d_fids[l_ids[i]] = d_ids[i];
2222
2223 }
2224
2225 /*---------------------------------------------------------------------------
2226 * Update mesh->i_face_cells array for ghost cells in standard halo.
2227 *
2228 * parameters:
2229 * mesh <-> pointer to cs_mesh_t structure
2230 * face_ifs <-- faces interface
2231 * cell_faces_idx <-- "cell -> faces" connectivity index
2232 * cell_faces_lst <-- "cell -> faces" connectivity list
2233 *---------------------------------------------------------------------------*/
2234
2235 static void
_update_i_face_cells(cs_mesh_t * mesh,cs_interface_set_t * face_ifs,cs_lnum_t * cell_faces_idx,cs_lnum_t * cell_faces_lst)2236 _update_i_face_cells(cs_mesh_t *mesh,
2237 cs_interface_set_t *face_ifs,
2238 cs_lnum_t *cell_faces_idx,
2239 cs_lnum_t *cell_faces_lst)
2240 {
2241 int rank, n_interfaces, distant_rank;
2242 cs_lnum_t i, j, gcell_id, face_id, if_id, shift, index_end;
2243 cs_lnum_t tr_id, start, end, length, mpi_count;
2244
2245 int local_rank_id = (cs_glob_n_ranks == 1) ? 0 : -1;
2246 cs_halo_t *halo = mesh->halo;
2247 cs_lnum_t *gcell_face_count = NULL, *l2d_fids = NULL;
2248 cs_lnum_t *send_buffer = NULL, *recv_buffer = NULL, *buffer = NULL;
2249 cs_lnum_t *send_shift = NULL, *recv_shift = NULL, *halo2ifs_rank = NULL;
2250
2251 const int n_init_perio = mesh->n_init_perio;
2252 const int n_ranks = cs_glob_n_ranks;
2253 const int local_rank = (cs_glob_rank_id == -1) ? 0 : cs_glob_rank_id;
2254 const int n_c_domains = halo->n_c_domains;
2255 const cs_lnum_t *s_index = halo->send_index;
2256 const cs_lnum_t *s_gcell_ids = halo->send_list;
2257 const cs_interface_t *interface = NULL;
2258 const cs_lnum_t *loc_ids = NULL;
2259 const cs_lnum_t *dist_ids = NULL;
2260 const cs_lnum_t *tr_index = NULL;
2261
2262 #if defined(HAVE_MPI)
2263 MPI_Request _request[128];
2264 MPI_Request *request = _request;
2265 MPI_Status _status[128];
2266 MPI_Status *status = _status;
2267
2268 if (n_c_domains*2 > 128) {
2269 BFT_MALLOC(request, n_c_domains*2, MPI_Request);
2270 BFT_MALLOC(status, n_c_domains*2, MPI_Status);
2271 }
2272 #endif
2273
2274 /* Get an interface set struct. on faces */
2275
2276 assert(face_ifs != NULL);
2277 n_interfaces = cs_interface_set_size(face_ifs);
2278
2279 /* Allocate auxiliary buffers */
2280
2281 BFT_MALLOC(l2d_fids, mesh->n_i_faces, cs_lnum_t);
2282 BFT_MALLOC(halo2ifs_rank, n_c_domains, cs_lnum_t);
2283
2284 /* If there is an extended neighborhood, n_c_domains can be different
2285 from n_interfaces (c_rank with only vertices on the interface) */
2286
2287 for (i = 0; i < n_c_domains; i++)
2288 halo2ifs_rank[i] = -1; /* No face interface between ranks */
2289
2290 index_end = 1; /* Only purely parralel faces are taking into account */
2291 index_end += 2*n_init_perio ; /* 2 transform. by initial periodicity */
2292
2293 /* Identify interface id and communicating rank */
2294
2295 for (if_id = 0; if_id < n_interfaces; if_id++) {
2296
2297 interface = cs_interface_set_get(face_ifs, if_id);
2298 distant_rank = cs_interface_rank(interface);
2299
2300 for (i = 0; i < n_c_domains; i++)
2301 if (halo->c_domain_rank[i] == distant_rank)
2302 break;
2303 assert(i < n_c_domains);
2304 halo2ifs_rank[i] = if_id;
2305
2306 } /* End of loop on interfaces */
2307
2308 /* First exchange number of faces linked to each ghost cells */
2309
2310 BFT_MALLOC(send_buffer, halo->n_send_elts[CS_HALO_EXTENDED], cs_lnum_t);
2311 BFT_MALLOC(send_shift, n_c_domains + 1, cs_lnum_t);
2312
2313 send_shift[0] = 0;
2314 for (i = 0; i < halo->n_send_elts[CS_HALO_EXTENDED]; i++)
2315 send_buffer[i] = 0;
2316
2317 /* Loop on communicating ranks to build send_buffer */
2318
2319 cs_interface_set_add_match_ids(face_ifs);
2320
2321 for (rank = 0, shift = 0; rank < n_c_domains; rank++) {
2322
2323 if_id = halo2ifs_rank[rank];
2324
2325 if (if_id > -1) { /* This communicating rank shares elements
2326 in the standard neighborhood */
2327
2328 interface = cs_interface_set_get(face_ifs, if_id);
2329 loc_ids = cs_interface_get_elt_ids(interface);
2330 dist_ids = cs_interface_get_match_ids(interface);
2331 tr_index = cs_interface_get_tr_index(interface);
2332
2333 assert(cs_interface_rank(interface) == halo->c_domain_rank[rank]);
2334
2335 /* Initalize l2d_fids */
2336
2337 for (i = 0; i < mesh->n_i_faces; i++)
2338 l2d_fids[i] = -1;
2339
2340 for (tr_id = 0; tr_id < index_end; tr_id++) {
2341
2342 _get_start_end_idx(mesh, s_index, halo->send_perio_lst,
2343 rank, tr_id, 0, /* STANDARD */
2344 &start, &end);
2345
2346 if (tr_id == 0) { /* Purely parallel faces */
2347
2348 if (tr_index != NULL) {
2349 _tag_interface_faces(-1, 0, tr_index[1],
2350 loc_ids, dist_ids, l2d_fids);
2351 assert(n_init_perio > 0);
2352 }
2353 else {
2354 _tag_interface_faces(-1, 0, cs_interface_size(interface),
2355 loc_ids, dist_ids, l2d_fids);
2356 assert(n_init_perio == 0);
2357 }
2358
2359 }
2360 else /* Periodic faces */
2361 _tag_interface_faces(tr_index[tr_id-1],
2362 tr_index[tr_id],
2363 tr_index[tr_id+1],
2364 loc_ids, dist_ids, l2d_fids);
2365
2366 for (i = start; i < end; i++) {
2367
2368 gcell_id = s_gcell_ids[i];
2369 for (j = cell_faces_idx[gcell_id];
2370 j < cell_faces_idx[gcell_id+1];
2371 j++) {
2372
2373 face_id = cell_faces_lst[j];
2374 if (face_id > -1) {
2375 if (l2d_fids[face_id] > -1) {
2376 shift++;
2377 send_buffer[i] += 1;
2378 }
2379 }
2380
2381 } /* End of loop cell -> face connectivity */
2382
2383 } /* End of loop on standard neighborhood for this tr_id and rank */
2384
2385 } /* End of loop on tr_ids */
2386
2387 } /* if_id > -1 */
2388
2389 send_shift[rank+1] = shift;
2390
2391 } /* End of loop on ranks */
2392
2393 /* Exchange data over the ranks => build face_count */
2394
2395 BFT_MALLOC(gcell_face_count, halo->n_elts[CS_HALO_EXTENDED], cs_lnum_t);
2396
2397 for (i = 0; i < halo->n_elts[CS_HALO_EXTENDED]; i++)
2398 gcell_face_count[i] = 0;
2399
2400 #if defined(HAVE_MPI)
2401 if (n_ranks > 1) {
2402
2403 /* Receive data from distant ranks */
2404
2405 mpi_count = 0;
2406
2407 for (rank = 0; rank < n_c_domains; rank++) {
2408
2409 start = halo->index[2*rank];
2410 length = halo->index[2*rank+1] - halo->index[2*rank];
2411
2412 if (halo->c_domain_rank[rank] != local_rank) {
2413
2414 buffer = gcell_face_count + start;
2415
2416 MPI_Irecv(buffer,
2417 length,
2418 CS_MPI_LNUM,
2419 halo->c_domain_rank[rank],
2420 halo->c_domain_rank[rank],
2421 cs_glob_mpi_comm,
2422 &(request[mpi_count++]));
2423
2424 }
2425 else
2426 local_rank_id = rank;
2427
2428 }
2429
2430 /* We wait for posting all receives (often recommended) */
2431
2432 MPI_Barrier(cs_glob_mpi_comm);
2433
2434 /* Send data to distant ranks */
2435
2436 for (rank = 0; rank < n_c_domains; rank++) {
2437
2438 /* If this is not the local rank */
2439
2440 if (halo->c_domain_rank[rank] != local_rank) {
2441
2442 start = s_index[2*rank];
2443 length = s_index[2*rank + 1] - s_index[2*rank];
2444 buffer = send_buffer + start;
2445
2446 MPI_Isend(buffer,
2447 length,
2448 CS_MPI_LNUM,
2449 halo->c_domain_rank[rank],
2450 local_rank,
2451 cs_glob_mpi_comm,
2452 &(request[mpi_count++]));
2453
2454 }
2455
2456 }
2457
2458 /* Wait for all exchanges */
2459
2460 MPI_Waitall(mpi_count, request, status);
2461 }
2462
2463 #endif /* defined(HAVE_MPI) */
2464
2465 /* Copy local values in case of periodicity */
2466
2467 if (n_init_perio > 0 && local_rank_id > -1) {
2468
2469 buffer = gcell_face_count + halo->index[2*local_rank_id];
2470 start = s_index[2*local_rank_id];
2471 length = s_index[2*local_rank_id + 1] - s_index[2*local_rank_id];
2472
2473 for (i = 0; i < length; i++)
2474 buffer[i] = send_buffer[start + i];
2475
2476 }
2477
2478 /* Build recv_shift */
2479
2480 BFT_MALLOC(recv_shift, n_c_domains + 1, cs_lnum_t);
2481
2482 recv_shift[0] = 0;
2483 for (rank = 0; rank < n_c_domains; rank++) {
2484 recv_shift[rank+1] = 0;
2485 for (i = halo->index[2*rank]; i < halo->index[2*rank+1]; i++)
2486 recv_shift[rank+1] += gcell_face_count[i];
2487 }
2488
2489 for (rank = 0; rank < n_c_domains; rank++)
2490 recv_shift[rank+1] += recv_shift[rank];
2491
2492 BFT_REALLOC(send_buffer, send_shift[n_c_domains], cs_lnum_t);
2493 BFT_MALLOC(recv_buffer, recv_shift[n_c_domains], cs_lnum_t);
2494
2495 /* Build send_buffer */
2496
2497 for (rank = 0; rank < n_c_domains; rank++) {
2498
2499 if_id = halo2ifs_rank[rank];
2500
2501 if (if_id > -1) {
2502
2503 interface = cs_interface_set_get(face_ifs, if_id);
2504 loc_ids = cs_interface_get_elt_ids(interface);
2505 dist_ids = cs_interface_get_match_ids(interface);
2506 tr_index = cs_interface_get_tr_index(interface);
2507 shift = send_shift[rank];
2508
2509 /* Initalize l2d_fids */
2510
2511 for (i = 0; i < mesh->n_i_faces; i++)
2512 l2d_fids[i] = -1;
2513
2514 for (tr_id = 0; tr_id < index_end; tr_id++) {
2515
2516 _get_start_end_idx(mesh, s_index, halo->send_perio_lst,
2517 rank, tr_id, 0, /* STANDARD */
2518 &start, &end);
2519
2520 if (tr_id == 0) { /* Purely parallel faces */
2521
2522 if (tr_index != NULL) {
2523 _tag_interface_faces(-1, 0, tr_index[1],
2524 loc_ids, dist_ids, l2d_fids);
2525 assert(n_init_perio > 0);
2526 }
2527 else {
2528 _tag_interface_faces(-1, 0, cs_interface_size(interface),
2529 loc_ids, dist_ids, l2d_fids);
2530 assert(n_init_perio == 0);
2531 }
2532
2533 }
2534 else /* Periodic faces */
2535 _tag_interface_faces(tr_index[tr_id-1],
2536 tr_index[tr_id],
2537 tr_index[tr_id+1],
2538 loc_ids, dist_ids, l2d_fids);
2539
2540 for (i = start; i < end; i++) {
2541
2542 gcell_id = s_gcell_ids[i];
2543
2544 for (j = cell_faces_idx[gcell_id];
2545 j < cell_faces_idx[gcell_id+1];
2546 j++) {
2547
2548 face_id = cell_faces_lst[j];
2549 if (face_id > -1)
2550 if (l2d_fids[face_id] > -1)
2551 send_buffer[shift++] = l2d_fids[face_id];
2552
2553 } /* End of loop cell -> face connectivity */
2554
2555 }
2556
2557 } /* End of loop on standard neighborhood for this communicating rank */
2558
2559 } /* if_id > -1 */
2560
2561 } /* End of loop on ranks */
2562
2563 /* Exchange face ids */
2564
2565 #if defined(HAVE_MPI)
2566 if (n_ranks > 1) { /* Exchange data over the ranks => build face_count */
2567
2568 /* Receive data from distant ranks */
2569
2570 mpi_count = 0;
2571
2572 for (rank = 0; rank < n_c_domains; rank++) {
2573
2574 start = recv_shift[rank];
2575 length = recv_shift[rank+1] - recv_shift[rank];
2576
2577 if (halo->c_domain_rank[rank] != local_rank) {
2578
2579 buffer = recv_buffer + start;
2580
2581 MPI_Irecv(buffer,
2582 length,
2583 CS_MPI_LNUM,
2584 halo->c_domain_rank[rank],
2585 halo->c_domain_rank[rank],
2586 cs_glob_mpi_comm,
2587 &(request[mpi_count++]));
2588
2589 }
2590
2591 }
2592
2593 /* We wait for posting all receives (often recommended) */
2594
2595 MPI_Barrier(cs_glob_mpi_comm);
2596
2597 /* Send data to distant ranks */
2598
2599 for (rank = 0; rank < n_c_domains; rank++) {
2600
2601 /* If this is not the local rank */
2602
2603 if (halo->c_domain_rank[rank] != local_rank) {
2604
2605 start = send_shift[rank];
2606 length = send_shift[rank + 1] - send_shift[rank];
2607 buffer = send_buffer + start;
2608
2609 MPI_Isend(buffer,
2610 length,
2611 CS_MPI_LNUM,
2612 halo->c_domain_rank[rank],
2613 local_rank,
2614 cs_glob_mpi_comm,
2615 &(request[mpi_count++]));
2616
2617 }
2618
2619 }
2620
2621 /* Wait for all exchanges */
2622
2623 MPI_Waitall(mpi_count, request, status);
2624 }
2625
2626 #endif /* defined(HAVE_MPI) */
2627
2628 /* Copy local values in case of periodicity */
2629
2630 if (n_init_perio > 0 && local_rank_id > -1) {
2631
2632 buffer = recv_buffer + recv_shift[local_rank_id];
2633 start = send_shift[local_rank_id];
2634 length = send_shift[local_rank_id + 1] - send_shift[local_rank_id];
2635
2636 for (i = 0; i < length; i++)
2637 buffer[i] = send_buffer[start + i];
2638
2639 }
2640
2641 /* Update face-> cells connect */
2642
2643 for (rank = 0; rank < n_c_domains; rank++) {
2644
2645 shift = recv_shift[rank];
2646
2647 for (i = halo->index[2*rank]; i < halo->index[2*rank+1]; i++) {
2648
2649 for (j = 0; j < gcell_face_count[i]; j++) {
2650
2651 face_id = recv_buffer[shift++];
2652
2653 /* In case of periodicity, a ghost cell may appear several times
2654 in the same halo. So we can have a face->cell connect up-to-date
2655 when it is not the first pass */
2656
2657 if (mesh->i_face_cells[face_id][0] == -1)
2658 mesh->i_face_cells[face_id][0] = mesh->n_cells + i;
2659 else if (mesh->i_face_cells[face_id][1] == -1)
2660 mesh->i_face_cells[face_id][1] = mesh->n_cells + i;
2661
2662 } /* End of loop on related faces */
2663
2664 } /* End of loop on standard neighborhood for this communicating rank */
2665
2666 } /* End of loop on ranks */
2667
2668 /* Free memory */
2669
2670 cs_interface_set_free_match_ids(face_ifs);
2671
2672 BFT_FREE(recv_buffer);
2673 BFT_FREE(send_buffer);
2674 BFT_FREE(send_shift);
2675 BFT_FREE(recv_shift);
2676 BFT_FREE(gcell_face_count);
2677 BFT_FREE(halo2ifs_rank);
2678 BFT_FREE(l2d_fids);
2679
2680 #if defined(HAVE_MPI)
2681 if (request != _request) {
2682 BFT_FREE(request);
2683 BFT_FREE(status);
2684 }
2685 #endif
2686
2687 /* Sanity check */
2688
2689 _check_i_face_cells(mesh);
2690 }
2691
2692 #if 0 /* TODO: check algorithm (deadlock on BG/L on one test case) */
2693
2694 /*---------------------------------------------------------------------------
2695 * Clean a halo i.e. delete rank(s) with no ghost cells.
2696 *
2697 * This case happens when an interface has only extended vertices
2698 * and the halo is standard.
2699 *
2700 * parameters:
2701 * mesh <-- pointer to a mesh structure
2702 *---------------------------------------------------------------------------*/
2703
2704 static void
2705 _clean_halo(cs_mesh_t *mesh)
2706 {
2707 cs_lnum_t rank_id, i, j;
2708
2709 cs_halo_t *halo = mesh->halo;
2710
2711 cs_lnum_t n_c_domains = halo->n_c_domains;
2712 cs_lnum_t n_real_c_domains = 0;
2713 cs_lnum_t counter = 0;
2714 cs_lnum_t *new_c_domain_rank = NULL;
2715 cs_lnum_t *new_perio_lst = NULL;
2716 cs_lnum_t *new_index = NULL;
2717
2718 const cs_lnum_t n_transforms = mesh->n_transforms;
2719
2720 cs_alloc_mode_t halo_alloc_mode = cs_halo_get_buffer_alloc_mode();
2721
2722 /* Is there something to do ? */
2723
2724 for (rank_id = 0; rank_id < n_c_domains; rank_id++)
2725 if (halo->send_index[2*rank_id+2] - halo->send_index[2*rank_id] > 0)
2726 n_real_c_domains++;
2727
2728 if (n_real_c_domains == n_c_domains)
2729 return;
2730
2731 /* halo->index, halo->perio_lst and n_c_domains need an update */
2732
2733 BFT_MALLOC(new_c_domain_rank, n_real_c_domains, cs_lnum_t);
2734
2735 CS_ALLOC_HD(new_index, 2*n_real_c_domains+1, cs_lnum_t, halo_alloc_mode);
2736
2737 if (n_transforms > 0)
2738 BFT_MALLOC(new_perio_lst, 4*n_transforms*n_real_c_domains, cs_lnum_t);
2739
2740 /* Define the new buffers */
2741
2742 new_index[0] = 0;
2743
2744 for (rank_id = 0; rank_id < n_c_domains; rank_id++) {
2745
2746 if (halo->send_index[2*rank_id+2] - halo->send_index[2*rank_id] > 0) {
2747
2748 new_c_domain_rank[counter] = halo->c_domain_rank[rank_id];
2749 new_index[2*counter+1] = halo->send_index[2*rank_id+1];
2750 new_index[2*counter+2] = halo->send_index[2*rank_id+2];
2751
2752 for (i = 0; i < n_transforms; i++)
2753 for (j = 0; j < 4; j++)
2754 new_perio_lst[4*counter + j + 4*n_real_c_domains*i]
2755 = halo->send_perio_lst[4*rank_id + j + 4*n_c_domains*i];
2756
2757 counter++;
2758
2759 } /* If there are elements for this rank */
2760
2761 } /* End of loop on frontier ranks */
2762
2763 /* Replace halo->send_perio_lst, halo->send_index and
2764 halo->c_domain_rank by new ones */
2765
2766 BFT_FREE(halo->c_domain_rank);
2767 CS_FREE_HD(halo->send_index);
2768
2769 if (n_transforms > 0)
2770 BFT_FREE(halo->send_perio_lst);
2771
2772 halo->n_c_domains = n_real_c_domains;
2773 halo->c_domain_rank = new_c_domain_rank;
2774
2775 halo->send_index = new_index;
2776 if (n_transforms > 0)
2777 halo->send_perio_lst = new_perio_lst;
2778
2779 /* Reallocation of halo's buffers */
2780
2781 BFT_REALLOC(halo->index, 2*n_real_c_domains+1, cs_lnum_t);
2782 BFT_REALLOC(halo->perio_lst, 4*n_transforms*n_real_c_domains, cs_lnum_t);
2783
2784 }
2785
2786 #endif /* #if 0 */
2787
2788 /*----------------------------------------------------------------------------
2789 * Define cell -> internal faces connectivity for ghost cells.
2790 *
2791 * Treatment of parallel and/or periodic halos for standard or extended
2792 * ghost cells according to halo type building option.
2793 *
2794 * parameters:
2795 * mesh <-- pointer to cs_mesh_t structure.
2796 * interface_set <-- pointer to cs_interface_set structure.
2797 * p_cell_faces_idx --> pointer to the connectivity index
2798 * p_cell_faces_lst --> pointer to the connectivity list
2799 *----------------------------------------------------------------------------*/
2800
2801 static void
_create_gcell_faces_connect(cs_mesh_t * mesh,const cs_interface_set_t * vertex_ifs,cs_lnum_t * p_cell_faces_idx[],cs_lnum_t * p_cell_faces_lst[])2802 _create_gcell_faces_connect(cs_mesh_t *mesh,
2803 const cs_interface_set_t *vertex_ifs,
2804 cs_lnum_t *p_cell_faces_idx[],
2805 cs_lnum_t *p_cell_faces_lst[])
2806 {
2807 cs_lnum_t i, fac_id, i_vtx, id1, id2, shift, vtx_id;
2808
2809 cs_lnum_t *vtx_tag = NULL;
2810 cs_lnum_t *cell_buffer = NULL, *cell_tag = NULL, *counter = NULL;
2811 cs_lnum_t *cell_faces_idx = NULL;
2812 cs_lnum_t *cell_faces_lst = NULL;
2813
2814 const cs_lnum_t n_i_faces = mesh->n_i_faces;
2815 const cs_lnum_t n_cells = mesh->n_cells;
2816 const cs_lnum_2_t *face_cells = (const cs_lnum_2_t *)(mesh->i_face_cells);
2817 const cs_lnum_t *fac_vtx_idx = mesh->i_face_vtx_idx;
2818 const cs_lnum_t *fac_vtx_lst = mesh->i_face_vtx_lst;
2819
2820 *p_cell_faces_idx = cell_faces_idx;
2821 *p_cell_faces_lst = cell_faces_lst;
2822
2823 if (vertex_ifs == NULL)
2824 return;
2825
2826 BFT_MALLOC(cell_faces_idx, n_cells+1, cs_lnum_t);
2827 BFT_MALLOC(cell_buffer, 2*n_cells, cs_lnum_t);
2828
2829 cell_tag = &(cell_buffer[0]);
2830 counter = &(cell_buffer[n_cells]);
2831
2832 cell_faces_idx[0] = 0;
2833 for (i = 0; i < n_cells; i++) {
2834 cell_faces_idx[i+1] = 0;
2835 cell_tag[i] = -1;
2836 }
2837
2838 assert(sizeof(cs_lnum_t) == sizeof(cs_lnum_t));
2839
2840 _get_vertex_tag(mesh->n_vertices, vertex_ifs, &vtx_tag);
2841
2842 for (fac_id = 0; fac_id < n_i_faces; fac_id++) {
2843
2844 for (i_vtx = fac_vtx_idx[fac_id];
2845 i_vtx < fac_vtx_idx[fac_id + 1];
2846 i_vtx++) {
2847
2848 vtx_id = fac_vtx_lst[i_vtx];
2849
2850 if (vtx_tag[vtx_id] == 1) {
2851
2852 id1 = face_cells[fac_id][0];
2853 id2 = face_cells[fac_id][1];
2854
2855 if (id1 < 0) {
2856 if (cell_tag[id2] != fac_id) {
2857 cell_tag[id2] = fac_id;
2858 cell_faces_idx[id2 + 1] += 1;
2859 }
2860 }
2861
2862 if (id2 < 0) {
2863 if (cell_tag[id1] != fac_id) {
2864 cell_tag[id1] = fac_id;
2865 cell_faces_idx[id1 + 1] += 1;
2866 }
2867 }
2868
2869 if (mesh->halo_type == CS_HALO_EXTENDED) {
2870 if (id1 >= 0 && id2 >= 0) {
2871
2872 if (cell_tag[id1] != fac_id) {
2873 cell_tag[id1] = fac_id;
2874 cell_faces_idx[id1 + 1] += 1;
2875 }
2876
2877 if (cell_tag[id2] != fac_id) {
2878 cell_tag[id2] = fac_id;
2879 cell_faces_idx[id2 + 1] += 1;
2880 }
2881
2882 }
2883 }
2884
2885 }
2886
2887 } /* End of loop on vertices */
2888
2889 } /* End of loop on internal faces */
2890
2891 /* Build index */
2892
2893 for (i = 0; i < n_cells; i++) {
2894 cell_faces_idx[i+1] += cell_faces_idx[i];
2895 counter[i] = 0;
2896 cell_tag[i] = -1;
2897 }
2898
2899 BFT_MALLOC(cell_faces_lst, cell_faces_idx[n_cells], cs_lnum_t);
2900
2901 for (fac_id = 0; fac_id < n_i_faces; fac_id++) {
2902
2903 for (i_vtx = fac_vtx_idx[fac_id];
2904 i_vtx < fac_vtx_idx[fac_id + 1];
2905 i_vtx++) {
2906
2907 vtx_id = fac_vtx_lst[i_vtx];
2908
2909 if (vtx_tag[vtx_id] == 1) {
2910
2911 id1 = face_cells[fac_id][0];
2912 id2 = face_cells[fac_id][1];
2913
2914 if (id1 < 0) {
2915 if (cell_tag[id2] != fac_id) {
2916
2917 cell_tag[id2] = fac_id;
2918 shift = cell_faces_idx[id2] + counter[id2];
2919 cell_faces_lst[shift] = fac_id;
2920 counter[id2] += 1;
2921
2922 }
2923 }
2924
2925 if (id2 < 0) {
2926 if (cell_tag[id1] != fac_id) {
2927
2928 cell_tag[id1] = fac_id;
2929 shift = cell_faces_idx[id1] + counter[id1];
2930 cell_faces_lst[shift] = fac_id;
2931 counter[id1] += 1;
2932
2933 }
2934 }
2935
2936 if (mesh->halo_type == CS_HALO_EXTENDED) {
2937 if (id1 >= 0 && id2 >= 0) {
2938
2939 if (cell_tag[id1] != fac_id) {
2940 cell_tag[id1] = fac_id;
2941 shift = cell_faces_idx[id1] + counter[id1];
2942 cell_faces_lst[shift] = fac_id;
2943 counter[id1] += 1;
2944 }
2945
2946 if (cell_tag[id2] != fac_id) {
2947 cell_tag[id2] = fac_id;
2948 shift = cell_faces_idx[id2] + counter[id2];
2949 cell_faces_lst[shift] = fac_id;
2950 counter[id2] += 1;
2951 }
2952
2953 }
2954 }
2955
2956 }
2957
2958 } /* End of loop on vertices */
2959
2960 } /* End of loop on internal faces */
2961
2962 BFT_FREE(vtx_tag);
2963 BFT_FREE(cell_buffer);
2964
2965 *p_cell_faces_idx = cell_faces_idx;
2966 *p_cell_faces_lst = cell_faces_lst;
2967 }
2968
2969 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
2970
2971 /*============================================================================
2972 * Public function definitions
2973 *============================================================================*/
2974
2975 /*----------------------------------------------------------------------------
2976 * Define halo structures for internal and distant ghost cells.
2977 *
2978 * parameters:
2979 * mesh <-- pointer to cs_mesh_t structure
2980 * face_ifs <-- pointer to faces interfaces
2981 * vertex_ifs <-- pointer to vertex interfaces
2982 * p_gcell_vtx_idx --> pointer to the connectivity index
2983 * p_gcell_vtx_lst --> pointer to the connectivity list
2984 *---------------------------------------------------------------------------*/
2985
2986 void
cs_mesh_halo_define(cs_mesh_t * mesh,cs_interface_set_t * face_ifs,cs_interface_set_t * vertex_ifs,cs_lnum_t * p_gcell_vtx_idx[],cs_lnum_t * p_gcell_vtx_lst[])2987 cs_mesh_halo_define(cs_mesh_t *mesh,
2988 cs_interface_set_t *face_ifs,
2989 cs_interface_set_t *vertex_ifs,
2990 cs_lnum_t *p_gcell_vtx_idx[],
2991 cs_lnum_t *p_gcell_vtx_lst[])
2992 {
2993 cs_lnum_t *send_gcell_vtx_idx = NULL, *send_gcell_vtx_lst = NULL;
2994 cs_lnum_t *gcell_vtx_idx = NULL, *gcell_vtx_lst = NULL;
2995 cs_lnum_t *gcell_faces_idx = NULL, *gcell_faces_lst = NULL;
2996 cs_halo_t *halo = mesh->halo;
2997
2998 halo->n_local_elts = mesh->n_cells;
2999
3000 /* Define cell -> internal faces connectivity for ghost cells */
3001
3002 _create_gcell_faces_connect(mesh,
3003 vertex_ifs,
3004 &gcell_faces_idx,
3005 &gcell_faces_lst);
3006
3007 /* Fill cs_halo_t structure for send_halo */
3008
3009 if (mesh->verbosity > 0) {
3010 bft_printf(_(" Local halo definition\n"));
3011 bft_printf_flush();
3012 }
3013
3014 _fill_send_halo(mesh,
3015 vertex_ifs,
3016 gcell_faces_idx,
3017 gcell_faces_lst);
3018
3019 #if 0
3020 /* Clean cs_halo_t structure.
3021 Remove communicating ranks with no ghost cells */
3022
3023 bft_printf(_(" Halo cleaning\n"));
3024 bft_printf_flush();
3025
3026 _clean_halo(mesh);
3027 #endif
3028
3029 /* Fill cs_halo_t structure for halo.
3030 We use the data from send_halo structure */
3031
3032 if (mesh->verbosity > 0) {
3033 bft_printf(_(" Distant halo creation\n"));
3034 bft_printf_flush();
3035 }
3036
3037 _fill_halo(mesh);
3038
3039 /* Update mesh structure elements bound to halo management */
3040
3041 mesh->n_ghost_cells = halo->n_elts[CS_HALO_EXTENDED];
3042 mesh->n_cells_with_ghosts = mesh->n_cells + mesh->n_ghost_cells;
3043
3044 /* Update connectivity between internal faces and cells */
3045
3046 if (cs_mesh_n_g_ghost_cells(mesh) > 0) {
3047
3048 /* Create local ghost connectivity for send_halo cells and send it.
3049 Receive ghost cells connectivity for halo cells. */
3050
3051 _create_send_gcell_vtx_connect(mesh,
3052 vertex_ifs,
3053 gcell_faces_idx,
3054 gcell_faces_lst,
3055 &send_gcell_vtx_idx,
3056 &send_gcell_vtx_lst);
3057
3058 _exchange_gcell_vtx_connect(mesh,
3059 send_gcell_vtx_idx,
3060 send_gcell_vtx_lst,
3061 &gcell_vtx_idx,
3062 &gcell_vtx_lst);
3063
3064 /* Free memory */
3065
3066 BFT_FREE(send_gcell_vtx_idx);
3067 BFT_FREE(send_gcell_vtx_lst);
3068
3069 /* Define mesh->i_face_cells array for ghost cells in standard halo and
3070 also ghost cells to ghost cells connectivity for standard and extended
3071 halo if necessary */
3072
3073 if (mesh->verbosity > 0) {
3074 bft_printf(_(" Updating the faces -> cells connectivity\n"));
3075 bft_printf_flush();
3076 }
3077
3078 _update_i_face_cells(mesh, face_ifs, gcell_faces_idx, gcell_faces_lst);
3079
3080 }
3081
3082 BFT_FREE(gcell_faces_idx);
3083 BFT_FREE(gcell_faces_lst);
3084
3085 *p_gcell_vtx_idx = gcell_vtx_idx;
3086 *p_gcell_vtx_lst = gcell_vtx_lst;
3087
3088 /* Update mesh structure elements bound to halo management */
3089
3090 if (mesh->n_ghost_cells > 0)
3091 BFT_REALLOC(mesh->cell_family, mesh->n_cells_with_ghosts, int);
3092
3093 #if 0 /* for debugging purposes */
3094 cs_halo_dump(halo, 1);
3095 #endif
3096 }
3097
3098 /*----------------------------------------------------------------------------*/
3099
3100 END_C_DECLS
3101