1 /*============================================================================
2 * Internal coupling: coupling for one instance of Code_Saturne
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 /*----------------------------------------------------------------------------
28 * Standard C library headers
29 *----------------------------------------------------------------------------*/
30
31 #include <assert.h>
32 #include <stdio.h>
33 #include <stdlib.h>
34 #include <string.h>
35 #include <limits.h>
36 #include <math.h>
37
38 #if defined(HAVE_MPI)
39 #include <mpi.h>
40 #endif
41
42 /*----------------------------------------------------------------------------
43 * Local headers
44 *----------------------------------------------------------------------------*/
45
46 #include "bft_mem.h"
47 #include "bft_printf.h"
48 #include "bft_error.h"
49
50 #include "fvm_defs.h"
51 #include "fvm_selector.h"
52
53 #include "cs_defs.h"
54 #include "cs_math.h"
55 #include "cs_sort.h"
56 #include "cs_search.h"
57 #include "cs_mesh_connect.h"
58 #include "cs_mesh_location.h"
59 #include "cs_coupling.h"
60 #include "cs_halo.h"
61 #include "cs_matrix.h"
62 #include "cs_mesh.h"
63 #include "cs_mesh_boundary.h"
64 #include "cs_mesh_group.h"
65 #include "cs_mesh_quantities.h"
66 #include "cs_convection_diffusion.h"
67 #include "cs_field.h"
68 #include "cs_field_operator.h"
69 #include "cs_selector.h"
70 #include "cs_parall.h"
71 #include "cs_prototypes.h"
72 #include "cs_velocity_pressure.h"
73
74 /*----------------------------------------------------------------------------
75 * Header for the current file
76 *----------------------------------------------------------------------------*/
77
78 #include "cs_internal_coupling.h"
79
80 /*----------------------------------------------------------------------------*/
81
82 BEGIN_C_DECLS
83
84 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
85
86 /*=============================================================================
87 * Local Macro Definitions
88 *============================================================================*/
89
90 /*============================================================================
91 * Local structure definitions
92 *============================================================================*/
93
94 /*============================================================================
95 * Static global variables
96 *============================================================================*/
97
98 static cs_internal_coupling_t *_internal_coupling = NULL;
99 static int _n_internal_couplings = 0;
100
101 /*============================================================================
102 * Prototypes for functions intended for use only by Fortran wrappers.
103 * (descriptions follow, with function bodies).
104 *============================================================================*/
105
106 void
107 cs_f_ic_field_coupled_faces(const int field_id,
108 bool **p);
109
110 /*============================================================================
111 * Fortran wrapper function definitions
112 *============================================================================*/
113
114 /*----------------------------------------------------------------------------
115 * Return a pointer to coupling faces indicator array by face id for a given
116 * field id.
117 *
118 * This function is intended for use by Fortran wrappers.
119 *
120 * parameters:
121 * field_id <-- field id
122 * p --> returned pointer
123 *----------------------------------------------------------------------------*/
124
125 void
cs_f_ic_field_coupled_faces(const int field_id,bool ** p)126 cs_f_ic_field_coupled_faces(const int field_id,
127 bool **p)
128 {
129 const cs_field_t* f = cs_field_by_id(field_id);
130 const int coupling_key_id = cs_field_key_id("coupling_entity");
131 int coupling_id = cs_field_get_key_int(f,
132 coupling_key_id);
133 const cs_internal_coupling_t *cpl
134 = cs_internal_coupling_by_id(coupling_id);
135
136 *p = cpl->coupled_faces;
137 }
138
139 /*=============================================================================
140 * Private function definitions
141 *============================================================================*/
142
143 /*----------------------------------------------------------------------------
144 * Return the equivalent heat transfer coefficient. If both terms are
145 * below a given tolerance, 0. is returned.
146 *
147 * parameters:
148 * h1 <-- first exchange coefficient
149 * h2 <-- second exchange coefficient
150 *
151 * return:
152 * value of equivalent exchange coefficient
153 *----------------------------------------------------------------------------*/
154
155 static inline cs_real_t
_calc_heq(cs_real_t h1,cs_real_t h2)156 _calc_heq(cs_real_t h1,
157 cs_real_t h2)
158 {
159 const cs_real_t h_eps = 1.e-12;
160
161 cs_real_t heq = 0.;
162 if (h1 + h2 > h_eps)
163 heq = h1 * h2 / (h1 + h2);
164
165 return heq;
166 }
167
168 /*----------------------------------------------------------------------------
169 * Compute the inverse of the face viscosity tensor and anisotropic vector
170 * taking into account the weight coefficients to update cocg for lsq gradient.
171 *
172 * parameters:
173 * wi <-- Weight coefficient of cell i
174 * wj <-- Weight coefficient of cell j
175 * d <-- IJ direction
176 * a <-- geometric weight J'F/I'J'
177 * ki_d --> Updated vector for cell i
178 *----------------------------------------------------------------------------*/
179
180 static inline void
_compute_ani_weighting_cocg(const cs_real_t wi[],const cs_real_t wj[],const cs_real_t d[],const cs_real_t a,cs_real_t ki_d[])181 _compute_ani_weighting_cocg(const cs_real_t wi[],
182 const cs_real_t wj[],
183 const cs_real_t d[],
184 const cs_real_t a,
185 cs_real_t ki_d[])
186 {
187 cs_real_t _d[3];
188 cs_real_6_t sum;
189 cs_real_6_t inv_wj;
190
191 for (int ii = 0; ii < 6; ii++)
192 sum[ii] = a*wi[ii] + (1. - a)*wj[ii];
193
194 cs_math_sym_33_inv_cramer(wj, inv_wj);
195
196 /* Note: K_i.K_f^-1 = SUM.K_j^-1
197 * K_j.K_f^-1 = SUM.K_i^-1
198 * So: K_i d = SUM.K_j^-1.IJ */
199
200 cs_math_sym_33_3_product(inv_wj, d, _d);
201 cs_math_sym_33_3_product(sum, _d, ki_d);
202 }
203
204 /*----------------------------------------------------------------------------
205 * Update R.H.S. for lsq gradient taking into account the weight coefficients.
206 *
207 * parameters:
208 * wi <-- Weight coefficient of cell i
209 * wj <-- Weight coefficient of cell j
210 * p_diff <-- R.H.S.
211 * d <-- R.H.S.
212 * a <-- geometric weight J'F/I'J'
213 * resi --> Updated R.H.S. for cell i
214 *----------------------------------------------------------------------------*/
215
216 static inline void
_compute_ani_weighting(const cs_real_t wi[],const cs_real_t wj[],const cs_real_t p_diff,const cs_real_t d[],const cs_real_t a,cs_real_t resi[])217 _compute_ani_weighting(const cs_real_t wi[],
218 const cs_real_t wj[],
219 const cs_real_t p_diff,
220 const cs_real_t d[],
221 const cs_real_t a,
222 cs_real_t resi[])
223 {
224 int ii;
225 cs_real_t _d[3];
226 cs_real_t ki_d[3];
227 cs_real_6_t inv_wj;
228 cs_real_6_t sum;
229
230 for (ii = 0; ii < 6; ii++)
231 sum[ii] = a*wi[ii] + (1. - a)*wj[ii];
232
233 cs_math_sym_33_inv_cramer(wj, inv_wj);
234
235 cs_math_sym_33_3_product(inv_wj, d, _d);
236 cs_math_sym_33_3_product(sum, _d, ki_d);
237
238 /* 1 / ||Ki. K_f^-1. IJ||^2 */
239 cs_real_t normi = 1. / cs_math_3_dot_product(ki_d, ki_d);
240
241 for (ii = 0; ii < 3; ii++)
242 resi[ii] += p_diff * ki_d[ii] * normi;
243 }
244
245 /*----------------------------------------------------------------------------
246 * Return the locator associated with given coupling entity and group number
247 *
248 * parameters:
249 * cpl <-- pointer to coupling structure
250 *----------------------------------------------------------------------------*/
251
252 static ple_locator_t*
_create_locator(cs_internal_coupling_t * cpl)253 _create_locator(cs_internal_coupling_t *cpl)
254 {
255 cs_lnum_t i, j;
256 cs_lnum_t ifac;
257
258 const cs_lnum_t n_local = cpl->n_local;
259 const int *c_tag = cpl->c_tag;
260 cs_lnum_t *faces_local_num = NULL;
261
262 fvm_nodal_t* nm = NULL;
263 int *tag_nm = NULL;
264 cs_lnum_t *faces_in_nm = NULL;
265
266 ple_coord_t* point_coords = NULL;
267
268 char mesh_name[16] = "locator";
269
270 /* Create locator */
271
272 #if defined(PLE_HAVE_MPI)
273 ple_locator_t *locator = ple_locator_create(cs_glob_mpi_comm,
274 cs_glob_n_ranks,
275 0);
276 #else
277 ple_locator_t *locator = ple_locator_create();
278 #endif
279
280 /* Create fvm_nodal_t structure */
281
282 BFT_MALLOC(faces_local_num, n_local, cs_lnum_t);
283 for (cs_lnum_t face_id = 0; face_id < n_local; face_id++)
284 faces_local_num[face_id] = cpl->faces_local[face_id] + 1;
285
286 nm = cs_mesh_connect_faces_to_nodal(cs_glob_mesh,
287 mesh_name,
288 false,
289 0,
290 n_local,
291 NULL,
292 faces_local_num); /* 1..n */
293
294 /* Tag fvm_nodal_t structure */
295
296 /* Number of faces to tag */
297 const int nfac_in_nm = fvm_nodal_get_n_entities(nm, 2);
298 /* Memory allocation */
299 BFT_MALLOC(faces_in_nm, nfac_in_nm, cs_lnum_t);
300 BFT_MALLOC(tag_nm, nfac_in_nm, int);
301 /* Get id of faces to tag in parent */
302 fvm_nodal_get_parent_num(nm, 2, faces_in_nm);
303 /* Tag faces */
304 for (cs_lnum_t ii = 0; ii < nfac_in_nm; ii++) {
305 /* Default tag is 0 */
306 tag_nm[ii] = 0;
307 for (cs_lnum_t jj = 0; jj < n_local; jj++) {
308 if (faces_in_nm[ii] == faces_local_num[jj]){
309 tag_nm[ii] = c_tag[jj];
310 break;
311 }
312 }
313 }
314 fvm_nodal_set_tag(nm, tag_nm, 2);
315 /* Free memory */
316 BFT_FREE(faces_in_nm);
317 BFT_FREE(tag_nm);
318 BFT_FREE(faces_local_num);
319
320 /* Creation of distant group cell centers */
321
322 BFT_MALLOC(point_coords, 3*n_local, cs_real_t);
323
324 for (i = 0; i < n_local; i++) {
325 ifac = cpl->faces_local[i]; /* 0..n-1 */
326 for (j = 0; j < 3; j++)
327 point_coords[3*i+j] = cs_glob_mesh_quantities->b_face_cog[3*ifac+j];
328 }
329
330 /* Locator initialization */
331
332 ple_locator_set_mesh(locator,
333 nm,
334 NULL,
335 0,
336 1.1, /* TODO */
337 cs_glob_mesh->dim,
338 n_local,
339 NULL,
340 c_tag,
341 point_coords,
342 NULL,
343 cs_coupling_mesh_extents,
344 cs_coupling_point_in_mesh_p);
345 /* Free memory */
346 nm = fvm_nodal_destroy(nm);
347 BFT_FREE(point_coords);
348 return locator;
349 }
350
351 /*----------------------------------------------------------------------------
352 * Destruction of given internal coupling structure.
353 *
354 * parameters:
355 * cpl <-> pointer to coupling structure to destroy
356 *----------------------------------------------------------------------------*/
357
358 static void
_destroy_entity(cs_internal_coupling_t * cpl)359 _destroy_entity(cs_internal_coupling_t *cpl)
360 {
361 BFT_FREE(cpl->c_tag);
362 BFT_FREE(cpl->faces_local);
363 BFT_FREE(cpl->faces_distant);
364 BFT_FREE(cpl->g_weight);
365 BFT_FREE(cpl->ci_cj_vect);
366 BFT_FREE(cpl->offset_vect);
367 BFT_FREE(cpl->coupled_faces);
368 BFT_FREE(cpl->cells_criteria);
369 BFT_FREE(cpl->faces_criteria);
370 BFT_FREE(cpl->interior_faces_group_name);
371 BFT_FREE(cpl->exterior_faces_group_name);
372 BFT_FREE(cpl->volume_zone_ids);
373 ple_locator_destroy(cpl->locator);
374 }
375
376 /*----------------------------------------------------------------------------
377 * Compute geometrical face weights around coupling interface.
378 *
379 * parameters:
380 * cpl <-- pointer to coupling structure
381 *----------------------------------------------------------------------------*/
382
383 static void
_compute_geometrical_face_weight(const cs_internal_coupling_t * cpl)384 _compute_geometrical_face_weight(const cs_internal_coupling_t *cpl)
385 {
386 cs_lnum_t face_id, cell_id;
387
388 const cs_lnum_t n_local = cpl->n_local;
389 const cs_lnum_t *faces_local = cpl->faces_local;
390 const cs_lnum_t n_distant = cpl->n_distant;
391 const cs_lnum_t *faces_distant = cpl->faces_distant;
392 const cs_real_3_t *ci_cj_vect = (const cs_real_3_t *)cpl->ci_cj_vect;
393 cs_real_t* g_weight = cpl->g_weight;
394
395 const cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
396 const cs_mesh_t *m = cs_glob_mesh;
397
398 const cs_real_t* cell_cen = fvq->cell_cen;
399 const cs_real_t* b_face_cog = fvq->b_face_cog;
400 const cs_real_t* diipb = fvq->diipb;
401 const cs_real_t* b_face_surf = cs_glob_mesh_quantities->b_face_surf;
402 const cs_real_t* b_face_normal = cs_glob_mesh_quantities->b_face_normal;
403
404 /* Store local FI' distances in gweight_distant */
405
406 cs_real_t *g_weight_distant = NULL;
407 BFT_MALLOC(g_weight_distant, n_distant, cs_real_t);
408 for (cs_lnum_t ii = 0; ii < n_distant; ii++) {
409 cs_real_t dv[3];
410 face_id = faces_distant[ii];
411 cell_id = m->b_face_cells[face_id];
412
413 for (cs_lnum_t jj = 0; jj < 3; jj++)
414 dv[jj] = - diipb[3*face_id + jj] - cell_cen[3*cell_id +jj]
415 + b_face_cog[3*face_id +jj];
416
417 g_weight_distant[ii] = cs_math_3_norm(dv);
418 }
419
420 /* Exchange FI' distances */
421
422 cs_internal_coupling_exchange_var(cpl,
423 1,
424 g_weight_distant,
425 g_weight);
426 /* Free memory */
427 BFT_FREE(g_weight_distant);
428
429 /* Normalise the distance to obtain geometrical weights */
430
431 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
432 face_id = faces_local[ii];
433 g_weight[ii] /= ( ci_cj_vect[ii][0]*b_face_normal[3*face_id+0]
434 + ci_cj_vect[ii][1]*b_face_normal[3*face_id+1]
435 + ci_cj_vect[ii][2]*b_face_normal[3*face_id+2])
436 / b_face_surf[face_id];
437 }
438 }
439
440 /*----------------------------------------------------------------------------
441 * Compute r_weight around coupling interface based on diffusivity c_weight.
442 *
443 * parameters:
444 * cpl <-- pointer to coupling structure
445 * c_weight <-- diffusivity
446 * r_weight --> physical face weight
447 *----------------------------------------------------------------------------*/
448
449 static void
_compute_physical_face_weight(const cs_internal_coupling_t * cpl,const cs_real_t c_weight[],cs_real_t rweight[])450 _compute_physical_face_weight(const cs_internal_coupling_t *cpl,
451 const cs_real_t c_weight[],
452 cs_real_t rweight[])
453 {
454 cs_lnum_t face_id, cell_id;
455
456 const cs_lnum_t n_local = cpl->n_local;
457 const cs_lnum_t *faces_local = cpl->faces_local;
458 const cs_real_t* g_weight = cpl->g_weight;
459
460 const cs_mesh_t* m = cs_glob_mesh;
461 const cs_lnum_t *restrict b_face_cells
462 = (const cs_lnum_t *restrict)m->b_face_cells;
463
464 /* Exchange c_weight */
465
466 cs_real_t *c_weight_local = NULL;
467 BFT_MALLOC(c_weight_local, n_local, cs_real_t);
468 cs_internal_coupling_exchange_by_cell_id(cpl,
469 1,
470 c_weight,
471 c_weight_local);
472
473 /* Compute rweight */
474
475 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
476 face_id = faces_local[ii];
477 cell_id = b_face_cells[face_id];
478 cs_real_t ki = c_weight[cell_id];
479 cs_real_t kj = c_weight_local[ii];
480 cs_real_t pond = g_weight[ii];
481 rweight[ii] = kj / ( pond * ki + (1. - pond) * kj);
482 }
483
484 /* Free memory */
485 BFT_FREE(c_weight_local);
486 }
487
488 /*----------------------------------------------------------------------------
489 * Compute offset vector on coupled faces
490 *
491 * parameters:
492 * cpl <-> pointer to coupling entity
493 *----------------------------------------------------------------------------*/
494
495 static void
_compute_offset(const cs_internal_coupling_t * cpl)496 _compute_offset(const cs_internal_coupling_t *cpl)
497 {
498 cs_lnum_t face_id, cell_id;
499
500 const cs_lnum_t n_local = cpl->n_local;
501 const cs_lnum_t *faces_local = cpl->faces_local;
502 const cs_real_t *g_weight = cpl->g_weight;
503 cs_real_3_t *offset_vect = cpl->offset_vect;
504
505 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
506 const cs_mesh_t *m = cs_glob_mesh;
507 const cs_real_t *cell_cen = mq->cell_cen;
508 const cs_real_t *b_face_cog = mq->b_face_cog;
509 const cs_lnum_t *b_face_cells = m->b_face_cells;
510
511 /* Exchange cell center location */
512
513 cs_real_t *cell_cen_local = NULL;
514 BFT_MALLOC(cell_cen_local, 3 * n_local, cs_real_t);
515 cs_internal_coupling_exchange_by_cell_id(cpl,
516 3,
517 mq->cell_cen,
518 cell_cen_local);
519
520 /* Compute OF vectors */
521
522 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
523 face_id = faces_local[ii];
524 cell_id = b_face_cells[face_id];
525
526 for (cs_lnum_t jj = 0; jj < 3; jj++) {
527 cs_real_t xxd = cell_cen_local[3*ii + jj];
528 cs_real_t xxl = cell_cen[3*cell_id + jj];
529 offset_vect[ii][jj] = b_face_cog[3*face_id + jj]
530 - ( g_weight[ii]*xxl
531 + (1. - g_weight[ii])*xxd);
532 }
533 }
534 /* Free memory */
535 BFT_FREE(cell_cen_local);
536 }
537
538 /*----------------------------------------------------------------------------
539 * Compute cell centers vectors IJ on coupled faces
540 *
541 * parameters:
542 * cpl <-- pointer to coupling entity
543 *----------------------------------------------------------------------------*/
544
545 static void
_compute_ci_cj_vect(const cs_internal_coupling_t * cpl)546 _compute_ci_cj_vect(const cs_internal_coupling_t *cpl)
547 {
548 cs_lnum_t face_id, cell_id;
549
550 const cs_lnum_t n_local = cpl->n_local;
551 const cs_lnum_t *faces_local = cpl->faces_local;
552 cs_real_3_t *ci_cj_vect = cpl->ci_cj_vect;
553
554 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
555 const cs_mesh_t *m = cs_glob_mesh;
556
557 const cs_real_t *cell_cen = mq->cell_cen;
558 const cs_lnum_t *b_face_cells = m->b_face_cells;
559
560 /* Exchange cell center location */
561
562 cs_real_t *cell_cen_local = NULL;
563 BFT_MALLOC(cell_cen_local, 3 * n_local, cs_real_t);
564 cs_internal_coupling_exchange_by_cell_id(cpl,
565 3, /* dimension */
566 mq->cell_cen,
567 cell_cen_local);
568
569 /* Compute IJ vectors */
570
571 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
572 face_id = faces_local[ii];
573 cell_id = b_face_cells[face_id];
574
575 for (cs_lnum_t jj = 0; jj < 3; jj++) {
576 cs_real_t xxd = cell_cen_local[3*ii + jj];
577 cs_real_t xxl = cell_cen[3*cell_id + jj];
578 ci_cj_vect[ii][jj] = xxd - xxl;
579 }
580 }
581
582 /* Free memory */
583 BFT_FREE(cell_cen_local);
584
585 /* Compute geometric weights and iterative reconstruction vector */
586
587 _compute_geometrical_face_weight(cpl);
588 _compute_offset(cpl);
589 }
590
591 /*----------------------------------------------------------------------------
592 * Define component coupled_faces[] of given coupling entity.
593 *
594 * parameters:
595 * cpl <-> pointer to coupling structure to modify
596 *----------------------------------------------------------------------------*/
597
598 static void
_initialize_coupled_faces(cs_internal_coupling_t * cpl)599 _initialize_coupled_faces(cs_internal_coupling_t *cpl)
600 {
601 const cs_lnum_t n_local = cpl->n_local;
602 const cs_lnum_t *faces_local = cpl->faces_local;
603 const cs_mesh_t *m = cs_glob_mesh;
604 bool *coupled_faces = cpl->coupled_faces;
605
606 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
607 coupled_faces[face_id] = false;
608 }
609
610 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
611 cs_lnum_t face_id = faces_local[ii];
612 coupled_faces[face_id] = true;
613 }
614 }
615
616 /*----------------------------------------------------------------------------
617 * Initialize to 0 or NULL most of the fields in given coupling structure.
618 *
619 * parameters:
620 * cpl --> pointer to coupling structure to initialize
621 *----------------------------------------------------------------------------*/
622
623 static void
_cpl_initialize(cs_internal_coupling_t * cpl)624 _cpl_initialize(cs_internal_coupling_t *cpl)
625 {
626 cpl->locator = NULL;
627 cpl->c_tag = NULL;
628 cpl->cells_criteria = NULL;
629 cpl->faces_criteria = NULL;
630 cpl->interior_faces_group_name = NULL;
631 cpl->exterior_faces_group_name = NULL;
632
633 cpl->n_volume_zones = 0;
634 cpl->volume_zone_ids = NULL;
635
636 cpl->n_local = 0;
637 cpl->faces_local = NULL; /* Coupling boundary faces, numbered 0..n-1 */
638
639 cpl->n_distant = 0;
640 cpl->faces_distant = NULL;
641
642 cpl->coupled_faces = NULL;
643
644 cpl->g_weight = NULL;
645 cpl->ci_cj_vect = NULL;
646 cpl->offset_vect = NULL;
647 }
648
649 /*----------------------------------------------------------------------------
650 * Initialize coupling criteria from strings.
651 *
652 * parameters:
653 * criteria_cells <-- string criteria for the first group of cells
654 * criteria_faces <-- string criteria for faces
655 * cpl --> pointer to coupling structure to initialize
656 *----------------------------------------------------------------------------*/
657
658 static void
_criteria_initialize(const char criteria_cells[],const char criteria_faces[],cs_internal_coupling_t * cpl)659 _criteria_initialize(const char criteria_cells[],
660 const char criteria_faces[],
661 cs_internal_coupling_t *cpl)
662 {
663 BFT_MALLOC(cpl->cells_criteria, strlen(criteria_cells)+1, char);
664 strcpy(cpl->cells_criteria, criteria_cells);
665
666 if (criteria_faces != NULL) {
667 BFT_MALLOC(cpl->faces_criteria, strlen(criteria_faces)+1, char);
668 strcpy(cpl->faces_criteria, criteria_faces);
669 }
670 }
671
672 /*----------------------------------------------------------------------------
673 * Define face to face mappings for internal couplings.
674 *
675 * The caller is responsible for freeing the list.
676 *
677 * parameters:
678 * cpl <-> pointer to internal coupling structure
679 * coupling_id <-- associated coupling id
680 *----------------------------------------------------------------------------*/
681
682 static void
_auto_group_name(cs_internal_coupling_t * cpl,int coupling_id)683 _auto_group_name(cs_internal_coupling_t *cpl,
684 int coupling_id)
685 {
686 char group_name[64];
687 snprintf(group_name, 63, "auto:internal_coupling_%d", coupling_id);
688 group_name[63] = '\0';
689 BFT_REALLOC(cpl->faces_criteria,
690 strlen(group_name)+1,
691 char);
692 strcpy(cpl->faces_criteria, group_name);
693 }
694
695 /*----------------------------------------------------------------------------
696 * Get selected cells list.
697 *
698 * parameters:
699 * m <-- pointer to mesh structure to modify
700 * cpl <-- pointer to coupling structure to modify
701 * n_cells --> number of associated cells
702 * cell_list --> associated cells list
703 *----------------------------------------------------------------------------*/
704
705 static void
_get_cell_list(cs_mesh_t * m,cs_internal_coupling_t * cpl,cs_lnum_t * n_cells,cs_lnum_t ** cell_list)706 _get_cell_list(cs_mesh_t *m,
707 cs_internal_coupling_t *cpl,
708 cs_lnum_t *n_cells,
709 cs_lnum_t **cell_list)
710 {
711 cs_lnum_t n_selected_cells = 0;
712 cs_lnum_t *selected_cells = NULL;
713
714 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
715
716 BFT_MALLOC(selected_cells, n_cells_ext, cs_lnum_t);
717
718 if (cpl->cells_criteria != NULL) {
719 cs_selector_get_cell_list(cpl->cells_criteria,
720 &n_selected_cells,
721 selected_cells);
722 }
723
724 /* For zone selection, zones may not be built yet, so use selection
725 mechanism directly. */
726
727 else if (cpl->n_volume_zones > 0) {
728
729 int *cell_flag;
730 BFT_MALLOC(cell_flag, n_cells_ext, int);
731 for (cs_lnum_t i = 0; i < n_cells_ext; i++)
732 cell_flag[i] = 0;
733
734 for (int i = 0; i < cpl->n_volume_zones; i++) {
735 const cs_zone_t *z = cs_volume_zone_by_id(cpl->volume_zone_ids[i]);
736 const char *criteria
737 = cs_mesh_location_get_selection_string(z->location_id);
738
739 if (criteria == NULL)
740 bft_error
741 (__FILE__, __LINE__, 0,
742 _("Only zones based on selection criteria strings "
743 "(not functions) are currently\n"
744 "supperted for the selection of internal coupling volumes.\n\n"
745 "This is not the case for zone: \"%s\"."), z->name);
746
747 cs_selector_get_cell_list(criteria, &n_selected_cells, selected_cells);
748
749 for (cs_lnum_t j = 0; j < n_selected_cells; j++)
750 cell_flag[selected_cells[j]] = 1;
751 }
752
753 n_selected_cells = 0;
754 for (cs_lnum_t i = 0; i < m->n_cells; i++) {
755 if (cell_flag[i] == 1) {
756 selected_cells[n_selected_cells] = i;
757 n_selected_cells++;
758 }
759 }
760
761 BFT_FREE(cell_flag);
762 }
763
764 BFT_REALLOC(selected_cells, n_selected_cells, cs_lnum_t);
765
766 *n_cells = n_selected_cells;
767 *cell_list = selected_cells;
768 }
769
770 /*----------------------------------------------------------------------------
771 * Initialize internal coupling and insert boundaries
772 * using cell selection criteria ONLY.
773 *
774 * parameters:
775 * m <-> pointer to mesh structure to modify
776 * cpl <-> pointer to coupling structure to modify
777 *----------------------------------------------------------------------------*/
778
779 static void
_volume_initialize_insert_boundary(cs_mesh_t * m,cs_internal_coupling_t * cpl)780 _volume_initialize_insert_boundary(cs_mesh_t *m,
781 cs_internal_coupling_t *cpl)
782 {
783 cs_lnum_t n_sel_cells = 0;
784 cs_lnum_t *sel_cells = NULL;
785
786 /* Selection of Volume zone using volumic selection criteria*/
787
788 _get_cell_list(m,
789 cpl,
790 &n_sel_cells,
791 &sel_cells);
792
793 int coupling_id = _n_internal_couplings - 1;
794
795 _auto_group_name(cpl, coupling_id);
796
797 cs_mesh_boundary_insert_separating_cells(m,
798 cpl->faces_criteria,
799 n_sel_cells,
800 sel_cells);
801
802 /* Select faces adjacent to volume zone and add appropriate group
803 so as to be able to easily extract separate sides */
804
805 {
806 cs_lnum_t n_sel_faces = 0;
807 cs_lnum_t *sel_faces_ext = NULL, *sel_faces_int = NULL;
808 int *cell_flag;
809
810 BFT_MALLOC(cell_flag, m->n_cells, int);
811 for (cs_lnum_t i = 0; i < m->n_cells; i++)
812 cell_flag[i] = 0;
813
814 for (cs_lnum_t i = 0; i < n_sel_cells; i++)
815 cell_flag[sel_cells[i]] = 1;
816
817 BFT_MALLOC(sel_faces_ext, m->n_b_faces, cs_lnum_t);
818 cs_selector_get_b_face_list(cpl->faces_criteria,
819 &n_sel_faces,
820 sel_faces_ext);
821
822 cs_lnum_t n_sel_int = 0, n_sel_ext = 0;
823 BFT_MALLOC(sel_faces_int, n_sel_faces, cs_lnum_t);
824
825 for (cs_lnum_t i = 0; i < n_sel_faces; i++) {
826 cs_lnum_t face_id = sel_faces_ext[i];
827 if (cell_flag[m->b_face_cells[face_id]]) {
828 sel_faces_ext[n_sel_ext] = face_id;
829 n_sel_ext++;
830 }
831 else {
832 sel_faces_int[n_sel_int] = face_id;
833 n_sel_int++;
834 }
835 }
836
837 BFT_FREE(cell_flag);
838
839 if (cpl->exterior_faces_group_name != NULL) {
840 cs_mesh_group_b_faces_add(m,
841 cpl->exterior_faces_group_name,
842 n_sel_ext,
843 sel_faces_ext);
844 }
845
846 if (cpl->interior_faces_group_name != NULL) {
847 cs_mesh_group_b_faces_add(m,
848 cpl->interior_faces_group_name,
849 n_sel_int,
850 sel_faces_int);
851 }
852
853 BFT_FREE(sel_faces_int);
854 BFT_FREE(sel_faces_ext);
855 }
856
857 BFT_FREE(sel_cells);
858 }
859
860 /*----------------------------------------------------------------------------
861 * Initialize internal coupling locators using cell and face selection criteria.
862 *
863 * parameters:
864 * m <-> pointer to mesh structure to modify
865 * cpl <-> pointer to coupling structure to modify
866 *----------------------------------------------------------------------------*/
867
868 static void
_volume_face_initialize(cs_mesh_t * m,cs_internal_coupling_t * cpl)869 _volume_face_initialize(cs_mesh_t *m,
870 cs_internal_coupling_t *cpl)
871 {
872 cs_lnum_t n_selected_cells = 0;
873 cs_lnum_t *selected_faces = NULL;
874 cs_lnum_t *selected_cells = NULL;
875
876 cs_lnum_t *cell_tag = NULL;
877
878 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
879
880 /* Selection of Volume zone using selection criteria */
881
882 _get_cell_list(m,
883 cpl,
884 &n_selected_cells,
885 &selected_cells);
886
887 /* Initialization */
888
889 BFT_MALLOC(cell_tag, n_cells_ext, cs_lnum_t);
890 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++)
891 cell_tag[cell_id] = 2;
892
893 /* Tag cells */
894
895 for (cs_lnum_t ii = 0; ii < n_selected_cells; ii++) {
896 cs_lnum_t cell_id = selected_cells[ii];
897 cell_tag[cell_id] = 1;
898 }
899 if (cs_glob_mesh->halo != NULL)
900 cs_halo_sync_num(cs_glob_mesh->halo, CS_HALO_STANDARD, cell_tag);
901
902 /* Free memory */
903
904 BFT_FREE(selected_cells);
905
906 /* Selection of the interface */
907
908 cs_lnum_t n_selected_faces = 0;
909
910 BFT_MALLOC(selected_faces, m->n_b_faces, cs_lnum_t);
911 cs_selector_get_b_face_list(cpl->faces_criteria,
912 &n_selected_faces,
913 selected_faces);
914
915 /* reorder selected faces */
916 {
917 cs_lnum_t n = 0;
918 char *b_face_flag;
919 BFT_MALLOC(b_face_flag, m->n_b_faces, char);
920 for (cs_lnum_t i = 0; i < m->n_b_faces; i++)
921 b_face_flag[i] = 0;
922 for (cs_lnum_t i = 0; i < n_selected_faces; i++)
923 b_face_flag[selected_faces[i]] = 1;
924 for (cs_lnum_t i = 0; i < m->n_b_faces; i++) {
925 if (b_face_flag[i] == 1)
926 selected_faces[n++] = i;
927 }
928 assert(n == n_selected_faces);
929 BFT_FREE(b_face_flag);
930 }
931
932 /* Prepare locator */
933
934 cpl->n_local = n_selected_faces; /* WARNING: only numerically
935 valid for conformal meshes */
936
937 BFT_MALLOC(cpl->faces_local, cpl->n_local, cs_lnum_t);
938 BFT_MALLOC(cpl->c_tag, cpl->n_local, int);
939
940 for (cs_lnum_t ii = 0; ii < cpl->n_local; ii++) {
941 cs_lnum_t face_id = selected_faces[ii];
942 cpl->faces_local[ii] = face_id;
943 cs_lnum_t cell_id = m->b_face_cells[face_id];
944 cpl->c_tag[ii] = cell_tag[cell_id];
945 }
946
947 /* Free memory */
948
949 BFT_FREE(selected_faces);
950 BFT_FREE(cell_tag);
951 }
952
953 /*----------------------------------------------------------------------------
954 * Initialize locators
955 *
956 * parameters:
957 * m <-> pointer to mesh structure to modify
958 * cpl <-> pointer to coupling structure to modify
959 *----------------------------------------------------------------------------*/
960
961 static void
_locator_initialize(cs_mesh_t * m,cs_internal_coupling_t * cpl)962 _locator_initialize(cs_mesh_t *m,
963 cs_internal_coupling_t *cpl)
964 {
965 /* Initialize locator */
966
967 cpl->locator = _create_locator(cpl);
968 cpl->n_distant = ple_locator_get_n_dist_points(cpl->locator);
969 BFT_MALLOC(cpl->faces_distant,
970 cpl->n_distant,
971 cs_lnum_t);
972 const cs_lnum_t *faces_distant_num
973 = ple_locator_get_dist_locations(cpl->locator);
974
975 /* From 1..n to 0..n-1 */
976 for (cs_lnum_t i = 0; i < cpl->n_distant; i++)
977 cpl->faces_distant[i] = faces_distant_num[i] - 1;
978
979 /* Geometric quantities */
980
981 BFT_MALLOC(cpl->g_weight, cpl->n_local, cs_real_t);
982 BFT_MALLOC(cpl->ci_cj_vect, cpl->n_local, cs_real_3_t);
983 BFT_MALLOC(cpl->offset_vect, cpl->n_local, cs_real_3_t);
984
985 _compute_ci_cj_vect(cpl);
986
987 BFT_MALLOC(cpl->coupled_faces, m->n_b_faces, bool);
988 }
989
990 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
991
992 /*============================================================================
993 * Public function definitions
994 *============================================================================*/
995
996 /*----------------------------------------------------------------------------*/
997 /*!
998 * \brief Return number of defined internal couplings.
999 *
1000 * \return number of internal couplings
1001 */
1002 /*----------------------------------------------------------------------------*/
1003
1004 int
cs_internal_coupling_n_couplings(void)1005 cs_internal_coupling_n_couplings(void)
1006 {
1007 return _n_internal_couplings;
1008 }
1009
1010 /*----------------------------------------------------------------------------*/
1011 /*!
1012 * \brief Define coupling volume using given selection criteria.
1013 *
1014 * Then, this volume must be separated from the rest of the domain with a wall.
1015 *
1016 * \param[in] criteria_cells criteria for the first group of cells
1017 * \param[in] criteria_faces criteria for faces to be joined
1018 */
1019 /*----------------------------------------------------------------------------*/
1020
1021 void
cs_internal_coupling_add(const char criteria_cells[],const char criteria_faces[])1022 cs_internal_coupling_add(const char criteria_cells[],
1023 const char criteria_faces[])
1024 {
1025 BFT_REALLOC(_internal_coupling,
1026 _n_internal_couplings + 1,
1027 cs_internal_coupling_t);
1028
1029 cs_internal_coupling_t *cpl = _internal_coupling + _n_internal_couplings;
1030
1031 cpl->id = _n_internal_couplings;
1032
1033 _cpl_initialize(cpl);
1034
1035 _criteria_initialize(criteria_cells, criteria_faces, cpl);
1036
1037 _n_internal_couplings++;
1038 }
1039
1040 /*----------------------------------------------------------------------------*/
1041 /*!
1042 * \brief Define coupling volume using given criteria. Then, this volume will
1043 * be separated from the rest of the domain with thin walls.
1044 *
1045 * \param[in] criteria_cells criteria for the first group of cells
1046 */
1047 /*----------------------------------------------------------------------------*/
1048
1049 void
cs_internal_coupling_add_volume(const char criteria_cells[])1050 cs_internal_coupling_add_volume(const char criteria_cells[])
1051 {
1052 if (_n_internal_couplings > 0)
1053 bft_error(__FILE__, __LINE__, 0,
1054 "Only one volume can be added in this version.");
1055
1056 BFT_REALLOC(_internal_coupling,
1057 _n_internal_couplings + 1,
1058 cs_internal_coupling_t);
1059
1060 cs_internal_coupling_t *cpl = _internal_coupling + _n_internal_couplings;
1061
1062 cpl->id = _n_internal_couplings;
1063
1064 _cpl_initialize(cpl);
1065
1066 _criteria_initialize(criteria_cells, NULL, cpl);
1067
1068 _n_internal_couplings++;
1069 }
1070
1071 /*----------------------------------------------------------------------------*/
1072 /*!
1073 * \brief Define coupling volume using given cs_zone_t. Then, this volume will
1074 * be separated from the rest of the domain with thin walls.
1075 *
1076 * \param[in] z pointer to cs_volume_zone_t
1077 */
1078 /*----------------------------------------------------------------------------*/
1079
1080 void
cs_internal_coupling_add_volume_zone(const cs_zone_t * z)1081 cs_internal_coupling_add_volume_zone(const cs_zone_t *z)
1082 {
1083 int z_ids[] = {z->id};
1084
1085 cs_internal_coupling_add_volume_zones(1, z_ids);
1086 }
1087
1088 /*----------------------------------------------------------------------------*/
1089 /*!
1090 * \brief Define coupling volume using given cs_zone_t. Then, this volume will
1091 * be separated from the rest of the domain with thin walls.
1092 *
1093 * \param[in] n_zones number of associated volume zones
1094 * \param[in] zone_ids ids of associated volume zones
1095 */
1096 /*----------------------------------------------------------------------------*/
1097
1098 void
cs_internal_coupling_add_volume_zones(int n_zones,const int zone_ids[])1099 cs_internal_coupling_add_volume_zones(int n_zones,
1100 const int zone_ids[])
1101 {
1102 if (_n_internal_couplings > 0)
1103 bft_error(__FILE__, __LINE__, 0,
1104 "Only one volume can be added in this version.");
1105
1106 BFT_REALLOC(_internal_coupling,
1107 _n_internal_couplings + 1,
1108 cs_internal_coupling_t);
1109
1110 cs_internal_coupling_t *cpl = _internal_coupling + _n_internal_couplings;
1111
1112 _cpl_initialize(cpl);
1113
1114 cpl->id = _n_internal_couplings;
1115
1116 cpl->n_volume_zones = n_zones;
1117 BFT_MALLOC(cpl->volume_zone_ids, n_zones, int);
1118
1119 for (int i = 0; i < n_zones; i++)
1120 cpl->volume_zone_ids[i] = zone_ids[i];
1121
1122 _n_internal_couplings++;
1123 }
1124
1125 /*----------------------------------------------------------------------------*/
1126 /*!
1127 * \brief Define internal coupling volume boundary group names.
1128 *
1129 * This is used only for internal couplings based on a separation of volumes
1130 * (cs_internal_coupling_add_volume, cs_internal_coupling_add_volume_zone,
1131 * cs_internal_coupling_add_volume_zones).
1132 *
1133 * The interior name is used for faces adjacent to the main volume, and
1134 * the exterio name for faces adjacent to the selected (exterior) volume.
1135 *
1136 * This allows filtering faces on each side of the boundary in a simpler manner.
1137 *
1138 * \param[in, out] cpl pointer to mesh structure to modify
1139 * \param[in] criteria_cells criteria for the first group of cells
1140 */
1141 /*----------------------------------------------------------------------------*/
1142
1143 void
cs_internal_coupling_add_boundary_groups(cs_internal_coupling_t * cpl,const char * interior_name,const char * exterior_name)1144 cs_internal_coupling_add_boundary_groups(cs_internal_coupling_t *cpl,
1145 const char *interior_name,
1146 const char *exterior_name)
1147 {
1148 if (cpl != NULL && interior_name != NULL) {
1149 BFT_REALLOC(cpl->interior_faces_group_name, strlen(interior_name) + 1, char);
1150 strcpy(cpl->interior_faces_group_name, interior_name);
1151 }
1152
1153 if (cpl != NULL && exterior_name != NULL) {
1154 BFT_REALLOC(cpl->exterior_faces_group_name, strlen(exterior_name) + 1, char);
1155 strcpy(cpl->exterior_faces_group_name, exterior_name);
1156 }
1157 }
1158
1159 /*----------------------------------------------------------------------------*/
1160 /*!
1161 * \brief Impose wall BCs to internal coupled faces if not yet defined.
1162 *
1163 * \param[in, out] bc_type face boundary condition type
1164 */
1165 /*----------------------------------------------------------------------------*/
1166
1167 void
cs_internal_coupling_bcs(int bc_type[])1168 cs_internal_coupling_bcs(int bc_type[])
1169 {
1170 cs_internal_coupling_t *cpl;
1171
1172 for (int cpl_id = 0; cpl_id < _n_internal_couplings; cpl_id++) {
1173 cpl = _internal_coupling + cpl_id;
1174
1175 const cs_lnum_t n_local = cpl->n_local;
1176 const cs_lnum_t *faces_local = cpl->faces_local;
1177
1178 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
1179 cs_lnum_t face_id = faces_local[ii];
1180 if (bc_type[face_id] == 0)
1181 bc_type[face_id] = CS_SMOOTHWALL;
1182 }
1183 }
1184 }
1185
1186 /*----------------------------------------------------------------------------*/
1187 /*!
1188 * \brief Add contribution from coupled faces (internal coupling) to
1189 * initialisation for iterative scalar gradient calculation.
1190 *
1191 * \param[in] cpl pointer to coupling entity
1192 * \param[in] c_weight weighted gradient coefficient variable, or NULL
1193 * \param[in] pvar variable
1194 * \param[in, out] grad gradient
1195 */
1196 /*----------------------------------------------------------------------------*/
1197
1198 void
cs_internal_coupling_initialize_scalar_gradient(const cs_internal_coupling_t * cpl,const cs_real_t c_weight[],const cs_real_t pvar[],cs_real_3_t * restrict grad)1199 cs_internal_coupling_initialize_scalar_gradient(
1200 const cs_internal_coupling_t *cpl,
1201 const cs_real_t c_weight[],
1202 const cs_real_t pvar[],
1203 cs_real_3_t *restrict grad)
1204 {
1205 cs_lnum_t face_id, cell_id;
1206
1207 const cs_lnum_t n_local = cpl->n_local;
1208 const cs_lnum_t *faces_local = cpl->faces_local;
1209 const cs_real_t* g_weight = cpl->g_weight;
1210 cs_real_t *r_weight = NULL;
1211
1212 const cs_mesh_t* m = cs_glob_mesh;
1213 const cs_lnum_t *restrict b_face_cells
1214 = (const cs_lnum_t *restrict)m->b_face_cells;
1215
1216 const cs_mesh_quantities_t* fvq = cs_glob_mesh_quantities;
1217 const cs_real_3_t *restrict b_f_face_normal
1218 = (const cs_real_3_t *restrict)fvq->b_f_face_normal;
1219
1220 /* Exchange pvar */
1221 cs_real_t *pvar_local = NULL;
1222 BFT_MALLOC(pvar_local, n_local, cs_real_t);
1223 cs_internal_coupling_exchange_by_cell_id(cpl,
1224 1,
1225 pvar,
1226 pvar_local);
1227
1228 /* Preliminary step in case of heterogenous diffusivity */
1229
1230 if (c_weight != NULL) {
1231 BFT_MALLOC(r_weight, n_local, cs_real_t);
1232 _compute_physical_face_weight(cpl,
1233 c_weight, /* diffusivity */
1234 r_weight); /* physical face weight */
1235 /* Redefinition of rweight :
1236 Before : (1-g_weight)*rweight <==> 1 - ktpond
1237 Modif : rweight = ktpond
1238 Scope of this modification is local */
1239 for (cs_lnum_t ii = 0; ii < n_local; ii++)
1240 r_weight[ii] = 1.0 - (1.0-g_weight[ii]) * r_weight[ii];
1241 }
1242
1243 /* Add contribution */
1244
1245 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
1246 face_id = faces_local[ii];
1247 cell_id = b_face_cells[face_id];
1248
1249 /* compared to _initialize_scalar_gradient :
1250 1 - rweight <==> 1 - ktpond
1251 g_weight <==> weight = alpha_ij
1252 pvar[cell_id] <==> pvar[ii]
1253 pvar_local[ii] <==> pvar[jj]
1254 b_f_face_normal <==> i_f_face_normal */
1255 cs_real_t pfaci = (c_weight == NULL) ?
1256 (1.0-g_weight[ii]) * (pvar_local[ii] - pvar[cell_id]) :
1257 (1.0-r_weight[ii]) * (pvar_local[ii] - pvar[cell_id]);
1258
1259 for (cs_lnum_t j = 0; j < 3; j++)
1260 grad[cell_id][j] += pfaci * b_f_face_normal[face_id][j];
1261
1262 }
1263 /* Free memory */
1264 if (c_weight != NULL) BFT_FREE(r_weight);
1265 BFT_FREE(pvar_local);
1266 }
1267
1268 /*----------------------------------------------------------------------------
1269 * Add contribution from coupled faces (internal coupling) to initialisation
1270 * for iterative vector gradient calculation
1271 *
1272 * parameters:
1273 * cpl <-- pointer to coupling entity
1274 * c_weight <-- weighted gradient coefficient variable, or NULL
1275 * pvar <-- variable
1276 * grad <-> gradient
1277 *----------------------------------------------------------------------------*/
1278
1279 void
cs_internal_coupling_initialize_vector_gradient(const cs_internal_coupling_t * cpl,const cs_real_t c_weight[],const cs_real_3_t pvar[],cs_real_33_t * restrict grad)1280 cs_internal_coupling_initialize_vector_gradient(
1281 const cs_internal_coupling_t *cpl,
1282 const cs_real_t c_weight[],
1283 const cs_real_3_t pvar[],
1284 cs_real_33_t *restrict grad)
1285 {
1286 cs_lnum_t face_id, cell_id;
1287
1288 const cs_lnum_t n_local = cpl->n_local;
1289 const cs_lnum_t *faces_local = cpl->faces_local;
1290 const cs_real_t* g_weight = cpl->g_weight;
1291 cs_real_t *r_weight = NULL;
1292
1293 const cs_mesh_t* m = cs_glob_mesh;
1294 const cs_lnum_t *restrict b_face_cells
1295 = (const cs_lnum_t *restrict)m->b_face_cells;
1296
1297 const cs_mesh_quantities_t* fvq = cs_glob_mesh_quantities;
1298 const cs_real_3_t *restrict b_f_face_normal
1299 = (const cs_real_3_t *restrict)fvq->b_f_face_normal;
1300
1301 /* Exchange pvar */
1302 cs_real_3_t *pvar_local = NULL;
1303 BFT_MALLOC(pvar_local, n_local, cs_real_3_t);
1304 cs_internal_coupling_exchange_by_cell_id(cpl,
1305 3,
1306 (const cs_real_t *)pvar,
1307 (cs_real_t *)pvar_local);
1308
1309 /* Preliminary step in case of heterogenous diffusivity */
1310
1311 if (c_weight != NULL) {
1312 BFT_MALLOC(r_weight, n_local, cs_real_t);
1313 _compute_physical_face_weight(cpl,
1314 c_weight, /* diffusivity */
1315 r_weight); /* physical face weight */
1316 /* Redefinition of rweight :
1317 Before : (1-g_weight)*rweight <==> 1 - ktpond
1318 Modif : rweight = ktpond
1319 Scope of this modification is local */
1320 for (cs_lnum_t ii = 0; ii < n_local; ii++)
1321 r_weight[ii] = 1.0 - (1.0-g_weight[ii]) * r_weight[ii];
1322 }
1323
1324 /* Add contribution */
1325
1326 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
1327 face_id = faces_local[ii];
1328 cell_id = b_face_cells[face_id];
1329
1330 /* compared to _initialize_scalar_gradient :
1331 1 - rweight <==> 1 - ktpond
1332 g_weight <==> weight = alpha_ij
1333 pvar[cell_id] <==> pvar[ii]
1334 pvar_local[ii] <==> pvar[jj]
1335 b_f_face_normal <==> i_f_face_normal */
1336
1337 for (cs_lnum_t i = 0; i < 3; i++) {
1338 cs_real_t pfaci = (c_weight == NULL) ?
1339 (1.0-g_weight[ii]) * (pvar_local[ii][i] - pvar[cell_id][i]) :
1340 (1.0-r_weight[ii]) * (pvar_local[ii][i] - pvar[cell_id][i]);
1341
1342 for (cs_lnum_t j = 0; j < 3; j++)
1343 grad[cell_id][i][j] += pfaci * b_f_face_normal[face_id][j];
1344
1345 }
1346 }
1347
1348 /* Free memory */
1349
1350 BFT_FREE(r_weight);
1351 BFT_FREE(pvar_local);
1352 }
1353
1354 /*----------------------------------------------------------------------------
1355 * Add contribution from coupled faces (internal coupling) to initialisation
1356 * for iterative symmetric tensor gradient calculation
1357 *
1358 * parameters:
1359 * cpl <-- pointer to coupling entity
1360 * c_weight <-- weighted gradient coefficient variable, or NULL
1361 * pvar <-- variable
1362 * grad <-> gradient
1363 *----------------------------------------------------------------------------*/
1364
1365 void
cs_internal_coupling_initialize_tensor_gradient(const cs_internal_coupling_t * cpl,const cs_real_t c_weight[],const cs_real_6_t pvar[],cs_real_63_t * restrict grad)1366 cs_internal_coupling_initialize_tensor_gradient(
1367 const cs_internal_coupling_t *cpl,
1368 const cs_real_t c_weight[],
1369 const cs_real_6_t pvar[],
1370 cs_real_63_t *restrict grad)
1371 {
1372 cs_lnum_t face_id, cell_id;
1373
1374 const cs_lnum_t n_local = cpl->n_local;
1375 const cs_lnum_t *faces_local = cpl->faces_local;
1376 const cs_real_t* g_weight = cpl->g_weight;
1377 cs_real_t *r_weight = NULL;
1378
1379 const cs_mesh_t* m = cs_glob_mesh;
1380 const cs_lnum_t *restrict b_face_cells
1381 = (const cs_lnum_t *restrict)m->b_face_cells;
1382
1383 const cs_mesh_quantities_t* fvq = cs_glob_mesh_quantities;
1384 const cs_real_3_t *restrict b_f_face_normal
1385 = (const cs_real_3_t *restrict)fvq->b_f_face_normal;
1386
1387 /* Exchange pvar */
1388 cs_real_6_t *pvar_local = NULL;
1389 BFT_MALLOC(pvar_local, n_local, cs_real_6_t);
1390 cs_internal_coupling_exchange_by_cell_id(cpl,
1391 6,
1392 (const cs_real_t *)pvar,
1393 (cs_real_t *)pvar_local);
1394
1395 /* Preliminary step in case of heterogenous diffusivity */
1396
1397 if (c_weight != NULL) {
1398 BFT_MALLOC(r_weight, n_local, cs_real_t);
1399 _compute_physical_face_weight(cpl,
1400 c_weight, /* diffusivity */
1401 r_weight); /* physical face weight */
1402 /* Redefinition of rweight :
1403 Before : (1-g_weight)*rweight <==> 1 - ktpond
1404 Modif : rweight = ktpond
1405 Scope of this modification is local */
1406 for (cs_lnum_t ii = 0; ii < n_local; ii++)
1407 r_weight[ii] = 1.0 - (1.0-g_weight[ii]) * r_weight[ii];
1408 }
1409
1410 /* Add contribution */
1411
1412 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
1413 face_id = faces_local[ii];
1414 cell_id = b_face_cells[face_id];
1415
1416 /* compared to _initialize_scalar_gradient :
1417 1 - rweight <==> 1 - ktpond
1418 g_weight <==> weight = alpha_ij
1419 pvar[cell_id] <==> pvar[ii]
1420 pvar_local[ii] <==> pvar[jj]
1421 b_f_face_normal <==> i_f_face_normal */
1422
1423 for (cs_lnum_t i = 0; i < 6; i++) {
1424 cs_real_t pfaci = (c_weight == NULL) ?
1425 (1.0-g_weight[ii]) * (pvar_local[ii][i] - pvar[cell_id][i]) :
1426 (1.0-r_weight[ii]) * (pvar_local[ii][i] - pvar[cell_id][i]);
1427
1428 for (cs_lnum_t j = 0; j < 3; j++)
1429 grad[cell_id][i][j] += pfaci * b_f_face_normal[face_id][j];
1430
1431 }
1432 }
1433
1434 /* Free memory */
1435
1436 BFT_FREE(r_weight);
1437 BFT_FREE(pvar_local);
1438 }
1439
1440 /*----------------------------------------------------------------------------*/
1441 /*!
1442 * \brief Add internal coupling rhs contribution for iterative gradient
1443 * calculation
1444 *
1445 * \param[in] cpl pointer to coupling entity
1446 * \param[in] c_weight weighted gradient coefficient variable, or NULL
1447 * \param[in] grad pointer to gradient
1448 * \param[in] pvar pointer to variable
1449 * \param[in, out] rhs pointer to rhs contribution
1450 */
1451 /*----------------------------------------------------------------------------*/
1452
1453 void
cs_internal_coupling_iterative_scalar_gradient(const cs_internal_coupling_t * cpl,const cs_real_t c_weight[],cs_real_3_t * restrict grad,const cs_real_t pvar[],cs_real_3_t rhs[])1454 cs_internal_coupling_iterative_scalar_gradient(
1455 const cs_internal_coupling_t *cpl,
1456 const cs_real_t c_weight[],
1457 cs_real_3_t *restrict grad,
1458 const cs_real_t pvar[],
1459 cs_real_3_t rhs[])
1460 {
1461 cs_lnum_t face_id, cell_id;
1462
1463 const cs_lnum_t n_local = cpl->n_local;
1464 const cs_lnum_t *faces_local = cpl->faces_local;
1465 const cs_real_3_t *offset_vect = (const cs_real_3_t *)cpl->offset_vect;
1466 const cs_real_t* g_weight = cpl->g_weight;
1467 cs_real_t *r_weight = NULL;
1468
1469 const cs_mesh_t* m = cs_glob_mesh;
1470 const cs_lnum_t *restrict b_face_cells
1471 = (const cs_lnum_t *restrict)m->b_face_cells;
1472
1473 const cs_mesh_quantities_t* fvq = cs_glob_mesh_quantities;
1474 const cs_real_3_t *restrict b_f_face_normal
1475 = (const cs_real_3_t *restrict)fvq->b_f_face_normal;
1476
1477 /* Exchange grad and pvar */
1478 cs_real_3_t *grad_local = NULL;
1479 BFT_MALLOC(grad_local, n_local, cs_real_3_t);
1480 cs_real_t *pvar_local = NULL;
1481 BFT_MALLOC(pvar_local, n_local, cs_real_t);
1482 cs_internal_coupling_exchange_by_cell_id(cpl,
1483 3,
1484 (const cs_real_t *)grad,
1485 (cs_real_t *)grad_local);
1486 cs_internal_coupling_exchange_by_cell_id(cpl,
1487 1,
1488 pvar,
1489 pvar_local);
1490
1491 /* Preliminary step in case of heterogenous diffusivity */
1492
1493 if (c_weight != NULL) { /* Heterogenous diffusivity */
1494 BFT_MALLOC(r_weight, n_local, cs_real_t);
1495 _compute_physical_face_weight(cpl,
1496 c_weight, /* diffusivity */
1497 r_weight); /* physical face weight */
1498 /* Redefinition of rweight_* :
1499 Before : (1-g_weight)*rweight <==> 1 - ktpond
1500 Modif : rweight = ktpond
1501 Scope of this modification is local */
1502 for (cs_lnum_t ii = 0; ii < n_local; ii++)
1503 r_weight[ii] = 1.0 - (1.0-g_weight[ii]) * r_weight[ii];
1504 }
1505
1506 /* Compute rhs */
1507
1508 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
1509 face_id = faces_local[ii];
1510 cell_id = b_face_cells[face_id];
1511
1512 /*
1513 Remark: \f$ \varia_\face = \alpha_\ij \varia_\celli
1514 + (1-\alpha_\ij) \varia_\cellj\f$
1515 but for the cell \f$ \celli \f$ we remove
1516 \f$ \varia_\celli \sum_\face \vect{S}_\face = \vect{0} \f$
1517 and for the cell \f$ \cellj \f$ we remove
1518 \f$ \varia_\cellj \sum_\face \vect{S}_\face = \vect{0} \f$
1519 */
1520
1521 /* Reconstruction part
1522 compared to _iterative_scalar_gradient :
1523 g_weight <==> weight = alpha_ij
1524 pvar[cell_id] <==> pvar[cell_id1]
1525 pvar_local[ii] <==> pvar[cell_id2]
1526 b_f_face_normal <==> i_f_face_normal */
1527 cs_real_t pfaci = 0.5;
1528 pfaci *= offset_vect[ii][0]*(grad_local[ii][0]+grad[cell_id][0])
1529 +offset_vect[ii][1]*(grad_local[ii][1]+grad[cell_id][1])
1530 +offset_vect[ii][2]*(grad_local[ii][2]+grad[cell_id][2]);
1531 if (c_weight != NULL)
1532 pfaci += (1.0-r_weight[ii]) * (pvar_local[ii] - pvar[cell_id]);
1533 else
1534 pfaci += (1.0-g_weight[ii]) * (pvar_local[ii] - pvar[cell_id]);
1535
1536 for (cs_lnum_t j = 0; j < 3; j++)
1537 rhs[cell_id][j] += pfaci * b_f_face_normal[face_id][j];
1538 }
1539
1540 BFT_FREE(r_weight);
1541 BFT_FREE(grad_local);
1542 BFT_FREE(pvar_local);
1543 }
1544
1545 /*----------------------------------------------------------------------------*/
1546 /*!
1547 * \brief Add internal coupling rhs contribution for iterative vector gradient
1548 * calculation
1549 *
1550 * \param[in] cpl pointer to coupling entity
1551 * \param[in] c_weight weighted gradient coefficient variable, or NULL
1552 * \param[in] grad pointer to gradient
1553 * \param[in] pvar pointer to variable
1554 * \param[in, out] rhs pointer to rhs contribution
1555 */
1556 /*----------------------------------------------------------------------------*/
1557
1558 void
cs_internal_coupling_iterative_vector_gradient(const cs_internal_coupling_t * cpl,const cs_real_t c_weight[],cs_real_33_t * restrict grad,const cs_real_3_t pvar[],cs_real_33_t rhs[])1559 cs_internal_coupling_iterative_vector_gradient(
1560 const cs_internal_coupling_t *cpl,
1561 const cs_real_t c_weight[],
1562 cs_real_33_t *restrict grad,
1563 const cs_real_3_t pvar[],
1564 cs_real_33_t rhs[])
1565 {
1566 cs_lnum_t face_id, cell_id;
1567
1568 const cs_lnum_t n_local = cpl->n_local;
1569 const cs_lnum_t *faces_local = cpl->faces_local;
1570 const cs_real_3_t *offset_vect = (const cs_real_3_t *)cpl->offset_vect;
1571 const cs_real_t* g_weight = cpl->g_weight;
1572 cs_real_t *r_weight = NULL;
1573
1574 const cs_mesh_t* m = cs_glob_mesh;
1575 const cs_lnum_t *restrict b_face_cells
1576 = (const cs_lnum_t *restrict)m->b_face_cells;
1577
1578 const cs_mesh_quantities_t* fvq = cs_glob_mesh_quantities;
1579 const cs_real_3_t *restrict b_f_face_normal
1580 = (const cs_real_3_t *restrict)fvq->b_f_face_normal;
1581
1582 /* Exchange grad and pvar */
1583 cs_real_33_t *grad_local = NULL;
1584 BFT_MALLOC(grad_local, n_local, cs_real_33_t);
1585 cs_real_3_t *pvar_local = NULL;
1586 BFT_MALLOC(pvar_local, n_local, cs_real_3_t);
1587 cs_internal_coupling_exchange_by_cell_id(cpl,
1588 9,
1589 (const cs_real_t *)grad,
1590 (cs_real_t *)grad_local);
1591 cs_internal_coupling_exchange_by_cell_id(cpl,
1592 3,
1593 (const cs_real_t *)pvar,
1594 (cs_real_t *)pvar_local);
1595
1596 /* Preliminary step in case of heterogenous diffusivity */
1597
1598 if (c_weight != NULL) { /* Heterogenous diffusivity */
1599 BFT_MALLOC(r_weight, n_local, cs_real_t);
1600 _compute_physical_face_weight(cpl,
1601 c_weight, /* diffusivity */
1602 r_weight); /* physical face weight */
1603 /* Redefinition of rweight_* :
1604 Before : (1-g_weight)*rweight <==> 1 - ktpond
1605 Modif : rweight = ktpond
1606 Scope of this modification is local */
1607 for (cs_lnum_t ii = 0; ii < n_local; ii++)
1608 r_weight[ii] = 1.0 - (1.0-g_weight[ii]) * r_weight[ii];
1609 }
1610
1611 /* Compute rhs */
1612
1613 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
1614 face_id = faces_local[ii];
1615 cell_id = b_face_cells[face_id];
1616
1617 /*
1618 Remark: \f$ \varia_\face = \alpha_\ij \varia_\celli
1619 + (1-\alpha_\ij) \varia_\cellj\f$
1620 but for the cell \f$ \celli \f$ we remove
1621 \f$ \varia_\celli \sum_\face \vect{S}_\face = \vect{0} \f$
1622 and for the cell \f$ \cellj \f$ we remove
1623 \f$ \varia_\cellj \sum_\face \vect{S}_\face = \vect{0} \f$
1624 */
1625
1626 /* Reconstruction part
1627 compared to _iterative_scalar_gradient:
1628 g_weight <==> weight = alpha_ij
1629 pvar[cell_id] <==> pvar[cell_id1]
1630 pvar_local[ii] <==> pvar[cell_id2]
1631 b_f_face_normal <==> i_f_face_normal */
1632
1633 for (cs_lnum_t i = 0; i < 3; i++) {
1634 cs_real_t pfaci = 0.5;
1635 pfaci *= offset_vect[ii][0]*(grad_local[ii][i][0]+grad[cell_id][i][0])
1636 +offset_vect[ii][1]*(grad_local[ii][i][1]+grad[cell_id][i][1])
1637 +offset_vect[ii][2]*(grad_local[ii][i][2]+grad[cell_id][i][2]);
1638 if (c_weight != NULL)
1639 pfaci += (1.0-r_weight[ii]) * (pvar_local[ii][i] - pvar[cell_id][i]);
1640 else
1641 pfaci += (1.0-g_weight[ii]) * (pvar_local[ii][i] - pvar[cell_id][i]);
1642
1643 for (cs_lnum_t j = 0; j < 3; j++)
1644 rhs[cell_id][i][j] += pfaci * b_f_face_normal[face_id][j];
1645
1646 }
1647 }
1648
1649 BFT_FREE(r_weight);
1650 BFT_FREE(grad_local);
1651 BFT_FREE(pvar_local);
1652 }
1653
1654 /*----------------------------------------------------------------------------*/
1655 /*!
1656 * \brief Add internal coupling rhs contribution for iterative tensor gradient
1657 * calculation
1658 *
1659 * \param[in] cpl pointer to coupling entity
1660 * \param[in] c_weight weighted gradient coefficient variable, or NULL
1661 * \param[in] grad pointer to gradient
1662 * \param[in] pvar pointer to variable
1663 * \param[in, out] rhs pointer to rhs contribution
1664 */
1665 /*----------------------------------------------------------------------------*/
1666
1667 void
cs_internal_coupling_iterative_tensor_gradient(const cs_internal_coupling_t * cpl,const cs_real_t c_weight[],cs_real_63_t * restrict grad,const cs_real_6_t pvar[],cs_real_63_t rhs[])1668 cs_internal_coupling_iterative_tensor_gradient(
1669 const cs_internal_coupling_t *cpl,
1670 const cs_real_t c_weight[],
1671 cs_real_63_t *restrict grad,
1672 const cs_real_6_t pvar[],
1673 cs_real_63_t rhs[])
1674 {
1675 cs_lnum_t face_id, cell_id;
1676
1677 const cs_lnum_t n_local = cpl->n_local;
1678 const cs_lnum_t *faces_local = cpl->faces_local;
1679 const cs_real_3_t *offset_vect = (const cs_real_3_t *)cpl->offset_vect;
1680 const cs_real_t* g_weight = cpl->g_weight;
1681 cs_real_t *r_weight = NULL;
1682
1683 const cs_mesh_t* m = cs_glob_mesh;
1684 const cs_lnum_t *restrict b_face_cells
1685 = (const cs_lnum_t *restrict)m->b_face_cells;
1686
1687 const cs_mesh_quantities_t* fvq = cs_glob_mesh_quantities;
1688 const cs_real_3_t *restrict b_f_face_normal
1689 = (const cs_real_3_t *restrict)fvq->b_f_face_normal;
1690
1691 /* Exchange grad and pvar */
1692 cs_real_63_t *grad_local = NULL;
1693 BFT_MALLOC(grad_local, n_local, cs_real_63_t);
1694 cs_real_6_t *pvar_local = NULL;
1695 BFT_MALLOC(pvar_local, n_local, cs_real_6_t);
1696 cs_internal_coupling_exchange_by_cell_id(cpl,
1697 18,
1698 (const cs_real_t *)grad,
1699 (cs_real_t *)grad_local);
1700 cs_internal_coupling_exchange_by_cell_id(cpl,
1701 6,
1702 (const cs_real_t *)pvar,
1703 (cs_real_t *)pvar_local);
1704
1705 /* Preliminary step in case of heterogenous diffusivity */
1706
1707 if (c_weight != NULL) { /* Heterogenous diffusivity */
1708 BFT_MALLOC(r_weight, n_local, cs_real_t);
1709 _compute_physical_face_weight(cpl,
1710 c_weight, /* diffusivity */
1711 r_weight); /* physical face weight */
1712 /* Redefinition of rweight_* :
1713 Before : (1-g_weight)*rweight <==> 1 - ktpond
1714 Modif : rweight = ktpond
1715 Scope of this modification is local */
1716 for (cs_lnum_t ii = 0; ii < n_local; ii++)
1717 r_weight[ii] = 1.0 - (1.0-g_weight[ii]) * r_weight[ii];
1718 }
1719
1720 /* Compute rhs */
1721
1722 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
1723 face_id = faces_local[ii];
1724 cell_id = b_face_cells[face_id];
1725
1726 /*
1727 Remark: \f$ \varia_\face = \alpha_\ij \varia_\celli
1728 + (1-\alpha_\ij) \varia_\cellj\f$
1729 but for the cell \f$ \celli \f$ we remove
1730 \f$ \varia_\celli \sum_\face \vect{S}_\face = \vect{0} \f$
1731 and for the cell \f$ \cellj \f$ we remove
1732 \f$ \varia_\cellj \sum_\face \vect{S}_\face = \vect{0} \f$
1733 */
1734
1735 /* Reconstruction part
1736 compared to _iterative_scalar_gradient :
1737 g_weight <==> weight = alpha_ij
1738 pvar[cell_id] <==> pvar[cell_id1]
1739 pvar_local[ii] <==> pvar[cell_id2]
1740 b_f_face_normal <==> i_f_face_normal */
1741 for (cs_lnum_t i = 0; i < 6; i++) {
1742 cs_real_t pfaci = 0.5;
1743 pfaci *= offset_vect[ii][0]*(grad_local[ii][i][0]+grad[cell_id][i][0])
1744 +offset_vect[ii][1]*(grad_local[ii][i][1]+grad[cell_id][i][1])
1745 +offset_vect[ii][2]*(grad_local[ii][i][2]+grad[cell_id][i][2]);
1746 if (c_weight != NULL)
1747 pfaci += (1.0-r_weight[ii]) * (pvar_local[ii][i] - pvar[cell_id][i]);
1748 else
1749 pfaci += (1.0-g_weight[ii]) * (pvar_local[ii][i] - pvar[cell_id][i]);
1750
1751 for (cs_lnum_t j = 0; j < 3; j++)
1752 rhs[cell_id][i][j] += pfaci * b_f_face_normal[face_id][j];
1753 }
1754 }
1755
1756 if (c_weight != NULL) BFT_FREE(r_weight);
1757 BFT_FREE(grad_local);
1758 BFT_FREE(pvar_local);
1759 }
1760
1761 /*----------------------------------------------------------------------------*/
1762 /*!
1763 * \brief Add internal coupling contribution for reconstruction of the
1764 * gradient of a scalar.
1765 *
1766 * \param[in] cpl pointer to coupling entity
1767 * \param[in] r_grad pointer to reconstruction gradient
1768 * \param[in, out] grad pointer to gradient to be reconstructed var
1769 */
1770 /*----------------------------------------------------------------------------*/
1771
1772 void
cs_internal_coupling_reconstruct_scalar_gradient(const cs_internal_coupling_t * cpl,cs_real_3_t r_grad[restrict],cs_real_3_t grad[])1773 cs_internal_coupling_reconstruct_scalar_gradient
1774 (const cs_internal_coupling_t *cpl,
1775 cs_real_3_t r_grad[restrict],
1776 cs_real_3_t grad[])
1777 {
1778 const cs_lnum_t n_local = cpl->n_local;
1779 const cs_lnum_t *faces_local = cpl->faces_local;
1780 const cs_real_3_t *offset_vect = (const cs_real_3_t *)cpl->offset_vect;
1781
1782 const cs_mesh_t* m = cs_glob_mesh;
1783 const cs_lnum_t *restrict b_face_cells
1784 = (const cs_lnum_t *restrict)m->b_face_cells;
1785
1786 const cs_mesh_quantities_t* fvq = cs_glob_mesh_quantities;
1787 const cs_real_3_t *restrict b_f_face_normal
1788 = (const cs_real_3_t *restrict)fvq->b_f_face_normal;
1789
1790 /* Exchange r_grad */
1791
1792 cs_real_3_t *r_grad_local = NULL;
1793 BFT_MALLOC(r_grad_local, n_local, cs_real_3_t);
1794 cs_internal_coupling_exchange_by_cell_id(cpl,
1795 3,
1796 (const cs_real_t *)r_grad,
1797 (cs_real_t *)r_grad_local);
1798
1799 /* Compute rhs */
1800
1801 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
1802 cs_lnum_t face_id = faces_local[ii];
1803 cs_lnum_t cell_id = b_face_cells[face_id];
1804
1805 /* Reconstruction part
1806 compared to _iterative_scalar_gradient :
1807 b_f_face_normal <==> i_f_face_normal */
1808 cs_real_t rfac = 0.5;
1809 rfac *= offset_vect[ii][0]*(r_grad_local[ii][0]+r_grad[cell_id][0])
1810 +offset_vect[ii][1]*(r_grad_local[ii][1]+r_grad[cell_id][1])
1811 +offset_vect[ii][2]*(r_grad_local[ii][2]+r_grad[cell_id][2]);
1812
1813 for (cs_lnum_t ll = 0; ll < 3; ll++)
1814 grad[cell_id][ll] += rfac * b_f_face_normal[face_id][ll];
1815 }
1816
1817 BFT_FREE(r_grad_local);
1818 }
1819
1820 /*----------------------------------------------------------------------------*/
1821 /*!
1822 * \brief Add internal coupling contribution for reconstruction of the
1823 * gradient of a vector.
1824 *
1825 * \param[in] cpl pointer to coupling entity
1826 * \param[in] r_grad pointer to reconstruction gradient
1827 * \param[in, out] grad pointer to gradient to be reconstructed var
1828 */
1829 /*----------------------------------------------------------------------------*/
1830
1831 void
cs_internal_coupling_reconstruct_vector_gradient(const cs_internal_coupling_t * cpl,cs_real_33_t * restrict r_grad,cs_real_33_t grad[])1832 cs_internal_coupling_reconstruct_vector_gradient(
1833 const cs_internal_coupling_t *cpl,
1834 cs_real_33_t *restrict r_grad,
1835 cs_real_33_t grad[])
1836 {
1837 const cs_lnum_t n_local = cpl->n_local;
1838 const cs_lnum_t *faces_local = cpl->faces_local;
1839 const cs_real_3_t *offset_vect = (const cs_real_3_t *)cpl->offset_vect;
1840
1841 const cs_mesh_t* m = cs_glob_mesh;
1842 const cs_lnum_t *restrict b_face_cells
1843 = (const cs_lnum_t *restrict)m->b_face_cells;
1844
1845 const cs_mesh_quantities_t* fvq = cs_glob_mesh_quantities;
1846 const cs_real_3_t *restrict b_f_face_normal
1847 = (const cs_real_3_t *restrict)fvq->b_f_face_normal;
1848
1849 /* Exchange r_grad */
1850
1851 cs_real_33_t *r_grad_local = NULL;
1852 BFT_MALLOC(r_grad_local, n_local, cs_real_33_t);
1853 cs_internal_coupling_exchange_by_cell_id(cpl,
1854 9,
1855 (const cs_real_t *)r_grad,
1856 (cs_real_t *)r_grad_local);
1857
1858 /* Compute rhs */
1859
1860 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
1861 cs_lnum_t face_id = faces_local[ii];
1862 cs_lnum_t cell_id = b_face_cells[face_id];
1863
1864 /* Reconstruction part
1865 compared to _iterative_scalar_gradient :
1866 b_f_face_normal <==> i_f_face_normal */
1867 for (cs_lnum_t i = 0; i < 3; i++) {
1868 cs_real_t rfac = 0.5;
1869 rfac *= offset_vect[ii][0]*(r_grad_local[ii][i][0]+r_grad[cell_id][i][0])
1870 +offset_vect[ii][1]*(r_grad_local[ii][i][1]+r_grad[cell_id][i][1])
1871 +offset_vect[ii][2]*(r_grad_local[ii][i][2]+r_grad[cell_id][i][2]);
1872
1873 for (cs_lnum_t j = 0; j < 3; j++)
1874 grad[cell_id][i][j] += rfac * b_f_face_normal[face_id][j];
1875 }
1876 }
1877
1878 BFT_FREE(r_grad_local);
1879 }
1880
1881 /*----------------------------------------------------------------------------*/
1882 /*!
1883 * \brief Add internal coupling contribution for reconstruction of the
1884 * gradient of a symmetric tensor.
1885 *
1886 * \param[in] cpl pointer to coupling entity
1887 * \param[in] r_grad pointer to reconstruction gradient
1888 * \param[in, out] grad pointer to gradient to be reconstructed var
1889 */
1890 /*----------------------------------------------------------------------------*/
1891
1892 void
cs_internal_coupling_reconstruct_tensor_gradient(const cs_internal_coupling_t * cpl,cs_real_63_t * restrict r_grad,cs_real_63_t grad[])1893 cs_internal_coupling_reconstruct_tensor_gradient(
1894 const cs_internal_coupling_t *cpl,
1895 cs_real_63_t *restrict r_grad,
1896 cs_real_63_t grad[])
1897 {
1898 const cs_lnum_t n_local = cpl->n_local;
1899 const cs_lnum_t *faces_local = cpl->faces_local;
1900 const cs_real_3_t *offset_vect = (const cs_real_3_t *)cpl->offset_vect;
1901
1902 const cs_mesh_t* m = cs_glob_mesh;
1903 const cs_lnum_t *restrict b_face_cells
1904 = (const cs_lnum_t *restrict)m->b_face_cells;
1905
1906 const cs_mesh_quantities_t* fvq = cs_glob_mesh_quantities;
1907 const cs_real_3_t *restrict b_f_face_normal
1908 = (const cs_real_3_t *restrict)fvq->b_f_face_normal;
1909
1910 /* Exchange r_grad */
1911
1912 cs_real_63_t *r_grad_local = NULL;
1913 BFT_MALLOC(r_grad_local, n_local, cs_real_63_t);
1914 cs_internal_coupling_exchange_by_cell_id(cpl,
1915 18,
1916 (const cs_real_t *)r_grad,
1917 (cs_real_t *)r_grad_local);
1918
1919 /* Compute rhs */
1920
1921 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
1922 cs_lnum_t face_id = faces_local[ii];
1923 cs_lnum_t cell_id = b_face_cells[face_id];
1924
1925 /* Reconstruction part
1926 compared to _iterative_scalar_gradient :
1927 b_f_face_normal <==> i_f_face_normal */
1928 for (cs_lnum_t i = 0; i < 6; i++) {
1929 cs_real_t rfac = 0.5;
1930 rfac *= offset_vect[ii][0]*(r_grad_local[ii][i][0]+r_grad[cell_id][i][0])
1931 +offset_vect[ii][1]*(r_grad_local[ii][i][1]+r_grad[cell_id][i][1])
1932 +offset_vect[ii][2]*(r_grad_local[ii][i][2]+r_grad[cell_id][i][2]);
1933
1934 for (cs_lnum_t j = 0; j < 3; j++)
1935 grad[cell_id][i][j] += rfac * b_f_face_normal[face_id][j];
1936 }
1937 }
1938
1939 BFT_FREE(r_grad_local);
1940 }
1941
1942 /*----------------------------------------------------------------------------*/
1943 /*!
1944 * \brief Add internal coupling rhs contribution for LSQ gradient calculation
1945 *
1946 * \param[in] cpl pointer to coupling entity
1947 * \param[in] c_weight weighted gradient coefficient variable, or NULL
1948 * \param[in] w_stride stride of weighting coefficient
1949 * \param[in, out] rhsv pointer to rhs contribution
1950 */
1951 /*----------------------------------------------------------------------------*/
1952
1953 void
cs_internal_coupling_lsq_scalar_gradient(const cs_internal_coupling_t * cpl,const cs_real_t c_weight[],const int w_stride,cs_real_4_t rhsv[])1954 cs_internal_coupling_lsq_scalar_gradient(
1955 const cs_internal_coupling_t *cpl,
1956 const cs_real_t c_weight[],
1957 const int w_stride,
1958 cs_real_4_t rhsv[])
1959 {
1960 cs_real_t pfac;
1961 cs_real_3_t dc, fctb;
1962
1963 const cs_lnum_t n_local = cpl->n_local;
1964 const cs_lnum_t *faces_local = cpl->faces_local;
1965 const cs_lnum_t n_distant = cpl->n_distant;
1966 const cs_lnum_t *faces_distant = cpl->faces_distant;
1967 const cs_real_t* g_weight = cpl->g_weight;
1968 const cs_real_3_t *ci_cj_vect = (const cs_real_3_t *)cpl->ci_cj_vect;
1969
1970 const cs_mesh_t* m = cs_glob_mesh;
1971 const cs_lnum_t *restrict b_face_cells
1972 = (const cs_lnum_t *restrict)m->b_face_cells;
1973
1974 /* Variables for cases w_stride = 1 or 6 */
1975 const bool scalar_diff = (c_weight != NULL && w_stride == 1);
1976 const bool tensor_diff = (c_weight != NULL && w_stride == 6);
1977 cs_real_t *weight = NULL;
1978
1979 /* Exchange pvar stored in rhsv[][3] */
1980
1981 cs_real_t *pvar_distant = NULL;
1982 BFT_MALLOC(pvar_distant, n_distant, cs_real_t);
1983 for (cs_lnum_t ii = 0; ii < n_distant; ii++) {
1984 cs_lnum_t face_id = faces_distant[ii];
1985 cs_lnum_t cell_id = b_face_cells[face_id];
1986 pvar_distant[ii] = rhsv[cell_id][3];
1987 }
1988 cs_real_t *pvar_local = NULL;
1989 BFT_MALLOC(pvar_local, n_local, cs_real_t);
1990 cs_internal_coupling_exchange_var(cpl,
1991 1,
1992 pvar_distant,
1993 pvar_local);
1994 /* Free memory */
1995 BFT_FREE(pvar_distant);
1996
1997 /* Preliminary step in case of heterogenous diffusivity */
1998
1999 if (c_weight != NULL) { /* Heterogenous diffusivity */
2000 if (tensor_diff) {
2001 BFT_MALLOC(weight, 6*n_local, cs_real_t);
2002 cs_internal_coupling_exchange_by_cell_id(cpl,
2003 6,
2004 c_weight,
2005 weight);
2006 } else {
2007 BFT_MALLOC(weight, n_local, cs_real_t);
2008 _compute_physical_face_weight(cpl,
2009 c_weight, /* diffusivity */
2010 weight); /* physical face weight */
2011 }
2012 }
2013
2014 /* Compute rhs */
2015
2016 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
2017 cs_lnum_t face_id = faces_local[ii];
2018 cs_lnum_t cell_id = b_face_cells[face_id];
2019 for (cs_lnum_t ll = 0; ll < 3; ll++)
2020 dc[ll] = ci_cj_vect[ii][ll];
2021
2022 if (tensor_diff) {
2023 /* (P_j - P_i)*/
2024 cs_real_t p_diff = (pvar_local[ii] - rhsv[cell_id][3]);
2025
2026 _compute_ani_weighting(&c_weight[6*cell_id],
2027 &weight[6*ii],
2028 p_diff,
2029 dc,
2030 g_weight[ii],
2031 &rhsv[cell_id][0]);
2032 } else if (scalar_diff) {
2033 /* (P_j - P_i) / ||d||^2 */
2034 pfac = (pvar_local[ii] - rhsv[cell_id][3])
2035 / (dc[0]*dc[0] + dc[1]*dc[1] + dc[2]*dc[2]);
2036
2037 for (cs_lnum_t ll = 0; ll < 3; ll++)
2038 fctb[ll] = dc[ll] * pfac;
2039
2040 /* Compared with _lsq_scalar_gradient, weight from
2041 * _compute_physical_face_weight already contains denom */
2042 for (cs_lnum_t ll = 0; ll < 3; ll++)
2043 rhsv[cell_id][ll] += weight[ii] * fctb[ll];
2044 } else {
2045 /* (P_j - P_i) / ||d||^2 */
2046 pfac = (pvar_local[ii] - rhsv[cell_id][3])
2047 / (dc[0]*dc[0] + dc[1]*dc[1] + dc[2]*dc[2]);
2048
2049 for (cs_lnum_t ll = 0; ll < 3; ll++)
2050 fctb[ll] = dc[ll] * pfac;
2051
2052 for (cs_lnum_t ll = 0; ll < 3; ll++)
2053 rhsv[cell_id][ll] += fctb[ll];
2054 }
2055
2056 }
2057
2058 /* Free memory */
2059 BFT_FREE(weight);
2060 BFT_FREE(pvar_local);
2061 }
2062
2063 /*----------------------------------------------------------------------------*/
2064 /*!
2065 * \brief Add internal coupling rhs contribution for LSQ gradient calculation
2066 *
2067 * \param[in] cpl pointer to coupling entity
2068 * \param[in] c_weight weighted gradient coefficient variable, or NULL
2069 * \param[in] w_stride stride of weighting coefficient
2070 * \param[in] pvar pointer to variable
2071 * \param[in, out] rhs pointer to rhs contribution
2072 */
2073 /*----------------------------------------------------------------------------*/
2074
2075 void
cs_internal_coupling_lsq_vector_gradient(const cs_internal_coupling_t * cpl,const cs_real_t c_weight[],const int w_stride,const cs_real_3_t pvar[],cs_real_33_t rhs[])2076 cs_internal_coupling_lsq_vector_gradient(
2077 const cs_internal_coupling_t *cpl,
2078 const cs_real_t c_weight[],
2079 const int w_stride,
2080 const cs_real_3_t pvar[],
2081 cs_real_33_t rhs[])
2082 {
2083 cs_real_t pfac;
2084 cs_real_3_t dc, fctb;
2085
2086 const cs_lnum_t n_local = cpl->n_local;
2087 const cs_lnum_t *faces_local = cpl->faces_local;
2088 const cs_lnum_t n_distant = cpl->n_distant;
2089 const cs_lnum_t *faces_distant = cpl->faces_distant;
2090 const cs_real_t* g_weight = cpl->g_weight;
2091 const cs_real_3_t *ci_cj_vect = (const cs_real_3_t *)cpl->ci_cj_vect;
2092
2093 const cs_mesh_t* m = cs_glob_mesh;
2094 const cs_lnum_t *restrict b_face_cells
2095 = (const cs_lnum_t *restrict)m->b_face_cells;
2096
2097 /* Variables for cases w_stride = 1 or 6 */
2098 const bool scalar_diff = (c_weight != NULL && w_stride == 1);
2099 const bool tensor_diff = (c_weight != NULL && w_stride == 6);
2100 cs_real_t *weight = NULL;
2101
2102 /* Exchange pvar */
2103
2104 cs_real_3_t *pvar_distant = NULL;
2105 BFT_MALLOC(pvar_distant, n_distant, cs_real_3_t);
2106 for (cs_lnum_t ii = 0; ii < n_distant; ii++) {
2107 cs_lnum_t face_id = faces_distant[ii];
2108 cs_lnum_t cell_id = b_face_cells[face_id];
2109 for (cs_lnum_t i = 0; i < 3; i++)
2110 pvar_distant[ii][i] = pvar[cell_id][i];
2111 }
2112 cs_real_3_t *pvar_local = NULL;
2113 BFT_MALLOC(pvar_local, n_local, cs_real_3_t);
2114 cs_internal_coupling_exchange_var(cpl,
2115 3,
2116 (cs_real_t *)pvar_distant,
2117 (cs_real_t *)pvar_local);
2118 /* Free memory */
2119 BFT_FREE(pvar_distant);
2120
2121 /* Preliminary step in case of heterogenous diffusivity */
2122
2123 if (c_weight != NULL) { /* Heterogenous diffusivity */
2124 if (tensor_diff) {
2125 BFT_MALLOC(weight, 6*n_local, cs_real_t);
2126 cs_internal_coupling_exchange_by_cell_id(cpl,
2127 6,
2128 c_weight,
2129 weight);
2130 } else {
2131 BFT_MALLOC(weight, n_local, cs_real_t);
2132 _compute_physical_face_weight(cpl,
2133 c_weight, /* diffusivity */
2134 weight); /* physical face weight */
2135 }
2136 }
2137
2138 /* Compute rhs */
2139
2140 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
2141 cs_lnum_t face_id = faces_local[ii];
2142 cs_lnum_t cell_id = b_face_cells[face_id];
2143 for (cs_lnum_t ll = 0; ll < 3; ll++)
2144 dc[ll] = ci_cj_vect[ii][ll];
2145
2146 if (tensor_diff) {//FIXME
2147 /* (P_j - P_i)*/
2148 for (cs_lnum_t i = 0; i < 3; i++) {
2149 cs_real_t p_diff = (pvar_local[ii][i] - pvar[cell_id][i]);
2150
2151 _compute_ani_weighting(&c_weight[6*cell_id],
2152 &weight[6*ii],
2153 p_diff,
2154 dc,
2155 g_weight[ii],
2156 rhs[cell_id][i]);
2157 }
2158 } else if (scalar_diff) {
2159
2160 for (cs_lnum_t i = 0; i < 3; i++) {
2161 /* (P_j - P_i) / ||d||^2 */
2162 pfac = (pvar_local[ii][i] - pvar[cell_id][i])
2163 / (dc[0]*dc[0] + dc[1]*dc[1] + dc[2]*dc[2]);
2164
2165 for (cs_lnum_t j = 0; j < 3; j++)
2166 fctb[j] = dc[j] * pfac;
2167
2168 /* Compared with _lsq_scalar_gradient, weight from
2169 * _compute_physical_face_weight already contains denom */
2170 for (cs_lnum_t j = 0; j < 3; j++)
2171 rhs[cell_id][i][j] += weight[ii] * fctb[j];
2172 }
2173 } else {
2174 for (cs_lnum_t i = 0; i < 3; i++) {
2175 /* (P_j - P_i) / ||d||^2 */
2176 pfac = (pvar_local[ii][i] - pvar[cell_id][i])
2177 / (dc[0]*dc[0] + dc[1]*dc[1] + dc[2]*dc[2]);
2178
2179 for (cs_lnum_t j = 0; j < 3; j++)
2180 fctb[j] = dc[j] * pfac;
2181
2182 /* Compared with _lsq_scalar_gradient, weight from
2183 * _compute_physical_face_weight already contains denom */
2184 for (cs_lnum_t j = 0; j < 3; j++)
2185 rhs[cell_id][i][j] += fctb[j];
2186 }
2187 }
2188
2189 }
2190 /* Free memory */
2191 if (c_weight != NULL) BFT_FREE(weight);
2192 BFT_FREE(pvar_local);
2193 }
2194
2195 /*----------------------------------------------------------------------------
2196 * Modify LSQ COCG matrix to include internal coupling
2197 *
2198 * parameters:
2199 * cpl <-- pointer to coupling entity
2200 * cocg <-> cocg matrix modified
2201 *----------------------------------------------------------------------------*/
2202
2203 void
cs_internal_coupling_lsq_cocg_contribution(const cs_internal_coupling_t * cpl,cs_real_6_t cocg[])2204 cs_internal_coupling_lsq_cocg_contribution(const cs_internal_coupling_t *cpl,
2205 cs_real_6_t cocg[])
2206 {
2207 const cs_lnum_t n_local = cpl->n_local;
2208 const cs_lnum_t *faces_local = cpl->faces_local;
2209 const cs_real_3_t *ci_cj_vect = (const cs_real_3_t *)cpl->ci_cj_vect;
2210
2211 const cs_mesh_t* m = cs_glob_mesh;
2212 const cs_lnum_t *restrict b_face_cells
2213 = (const cs_lnum_t *restrict)m->b_face_cells;
2214
2215 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
2216 cs_lnum_t face_id = faces_local[ii];
2217 cs_lnum_t cell_id = b_face_cells[face_id];
2218
2219 cs_real_t dc[3];
2220 for (cs_lnum_t ll = 0; ll < 3; ll++)
2221 dc[ll] = ci_cj_vect[ii][ll];
2222 cs_real_t ddc = 1. / (dc[0]*dc[0] + dc[1]*dc[1] + dc[2]*dc[2]);
2223
2224 cocg[cell_id][0] += dc[0]*dc[0]*ddc;
2225 cocg[cell_id][1] += dc[1]*dc[1]*ddc;
2226 cocg[cell_id][2] += dc[2]*dc[2]*ddc;
2227 cocg[cell_id][3] += dc[0]*dc[1]*ddc;
2228 cocg[cell_id][4] += dc[1]*dc[2]*ddc;
2229 cocg[cell_id][5] += dc[0]*dc[2]*ddc;
2230 }
2231 }
2232
2233 /*----------------------------------------------------------------------------
2234 * Modify LSQ COCG matrix to include internal coupling
2235 * when diffusivity is a tensor
2236 *
2237 * parameters:
2238 * cpl <-- pointer to coupling entity
2239 * c_weight <-- weigthing coefficients
2240 * cocg <-> cocg matrix modified
2241 *----------------------------------------------------------------------------*/
2242
2243 void
cs_internal_coupling_lsq_cocg_weighted(const cs_internal_coupling_t * cpl,const cs_real_t * c_weight,cs_real_6_t cocg[])2244 cs_internal_coupling_lsq_cocg_weighted(const cs_internal_coupling_t *cpl,
2245 const cs_real_t *c_weight,
2246 cs_real_6_t cocg[])
2247 {
2248 const cs_lnum_t n_local = cpl->n_local;
2249 const cs_lnum_t *faces_local = cpl->faces_local;
2250 const cs_real_t* g_weight = cpl->g_weight;
2251 const cs_real_3_t *ci_cj_vect = (const cs_real_3_t *)cpl->ci_cj_vect;
2252
2253 const cs_mesh_t* m = cs_glob_mesh;
2254 const cs_lnum_t *restrict b_face_cells
2255 = (const cs_lnum_t *restrict)m->b_face_cells;
2256
2257 /* Exchange c_weight */
2258 cs_real_t *cwgt_local = NULL;
2259 BFT_MALLOC(cwgt_local, 6*n_local, cs_real_t);
2260 cs_internal_coupling_exchange_by_cell_id(cpl,
2261 6,
2262 c_weight,
2263 cwgt_local);
2264
2265 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
2266 cs_lnum_t face_id = faces_local[ii];
2267 cs_lnum_t cell_id = b_face_cells[face_id];
2268 cs_real_t dc[3];
2269 for (cs_lnum_t ll = 0; ll < 3; ll++)
2270 dc[ll] = ci_cj_vect[ii][ll];
2271
2272 /* Reproduce _compute_ani_weighting_cocg */
2273 cs_real_t pond = g_weight[ii];
2274 cs_real_t dc_i[3] = {0., 0., 0.};
2275 _compute_ani_weighting_cocg(&c_weight[cell_id*6],
2276 &cwgt_local[ii*6],
2277 dc,
2278 pond,
2279 dc_i);
2280
2281 cs_real_t i_dci = 1./ (dc_i[0]*dc_i[0] + dc_i[1]*dc_i[1] + dc_i[2]*dc_i[2]);
2282
2283 cocg[cell_id][0] += dc_i[0] * dc_i[0] * i_dci;
2284 cocg[cell_id][1] += dc_i[1] * dc_i[1] * i_dci;
2285 cocg[cell_id][2] += dc_i[2] * dc_i[2] * i_dci;
2286 cocg[cell_id][3] += dc_i[0] * dc_i[1] * i_dci;
2287 cocg[cell_id][4] += dc_i[1] * dc_i[2] * i_dci;
2288 cocg[cell_id][5] += dc_i[0] * dc_i[2] * i_dci;
2289 }
2290
2291 BFT_FREE(cwgt_local);
2292 }
2293
2294 /*----------------------------------------------------------------------------
2295 * Modify iterative COCG matrix to include internal coupling
2296 *
2297 * parameters:
2298 * cpl <-- pointer to coupling entity
2299 * cocg <-> cocg matrix modified
2300 *----------------------------------------------------------------------------*/
2301
2302 void
cs_internal_coupling_it_cocg_contribution(const cs_internal_coupling_t * cpl,cs_real_33_t cocg[])2303 cs_internal_coupling_it_cocg_contribution(const cs_internal_coupling_t *cpl,
2304 cs_real_33_t cocg[])
2305 {
2306 const cs_lnum_t n_local = cpl->n_local;
2307 const cs_lnum_t *faces_local = cpl->faces_local;
2308 const cs_real_3_t *offset_vect = (const cs_real_3_t *)cpl->offset_vect;
2309
2310 const cs_mesh_t* m = cs_glob_mesh;
2311 /* const int n_cells_ext = m->n_cells_with_ghosts; */
2312 const cs_lnum_t *restrict b_face_cells
2313 = (const cs_lnum_t *restrict)m->b_face_cells;
2314 const cs_mesh_quantities_t* fvq = cs_glob_mesh_quantities;
2315 const cs_real_3_t *restrict b_f_face_normal
2316 = (const cs_real_3_t *restrict)fvq->b_f_face_normal;
2317 const cs_real_t *restrict cell_vol = fvq->cell_vol;
2318
2319 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
2320 cs_lnum_t face_id = faces_local[ii];
2321 cs_lnum_t cell_id = b_face_cells[face_id];
2322
2323 for (cs_lnum_t ll = 0; ll < 3; ll++) {
2324 for (cs_lnum_t mm = 0; mm < 3; mm++)
2325 cocg[cell_id][ll][mm] -= 0.5 * offset_vect[ii][ll]
2326 * b_f_face_normal[face_id][mm] / cell_vol[cell_id];
2327 }
2328 }
2329 }
2330
2331 /*----------------------------------------------------------------------------*/
2332 /*!
2333 * \brief Destruction of all internal coupling related structures.
2334 */
2335 /*----------------------------------------------------------------------------*/
2336
2337 void
cs_internal_coupling_finalize(void)2338 cs_internal_coupling_finalize(void)
2339 {
2340 cs_internal_coupling_t* cpl;
2341 for (int cpl_id = 0; cpl_id < _n_internal_couplings; cpl_id++) {
2342 cpl = _internal_coupling + cpl_id;
2343 _destroy_entity(cpl);
2344 }
2345 BFT_FREE(_internal_coupling);
2346 _n_internal_couplings = 0;
2347 }
2348
2349 /*----------------------------------------------------------------------------*/
2350 /*!
2351 * \brief Return the coupling associated with a given coupling_id.
2352 *
2353 * \param[in] coupling_id associated with a coupling entity
2354 *
2355 * \return pointer to associated coupling structure
2356 */
2357 /*----------------------------------------------------------------------------*/
2358
2359 cs_internal_coupling_t *
cs_internal_coupling_by_id(int coupling_id)2360 cs_internal_coupling_by_id(int coupling_id)
2361 {
2362 if (coupling_id > -1 && coupling_id < _n_internal_couplings)
2363 return _internal_coupling + coupling_id;
2364 else
2365 bft_error(__FILE__, __LINE__, 0,
2366 "coupling_id = %d provided is invalid", coupling_id);
2367 return (cs_internal_coupling_t*)NULL;
2368 }
2369
2370 /*----------------------------------------------------------------------------*/
2371 /*!
2372 * \brief Exchange quantities from distant to local
2373 * (update local using distant).
2374 *
2375 * \param[in] cpl pointer to coupling entity
2376 * \param[in] stride stride (e.g. 1 for double, 3 for interleaved coordinates)
2377 * \param[in] distant distant values, size coupling->n_distant
2378 * \param[out] local local values, size coupling->n_local
2379 */
2380 /*----------------------------------------------------------------------------*/
2381
2382 void
cs_internal_coupling_exchange_var(const cs_internal_coupling_t * cpl,int stride,cs_real_t distant[],cs_real_t local[])2383 cs_internal_coupling_exchange_var(const cs_internal_coupling_t *cpl,
2384 int stride,
2385 cs_real_t distant[],
2386 cs_real_t local[])
2387 {
2388 ple_locator_exchange_point_var(cpl->locator,
2389 distant,
2390 local,
2391 NULL,
2392 sizeof(cs_real_t),
2393 stride,
2394 0);
2395 }
2396
2397 /*----------------------------------------------------------------------------*/
2398 /*!
2399 * \brief Exchange variable between groups using cell id.
2400 *
2401 * \param[in] cpl pointer to coupling entity
2402 * \param[in] stride number of values (non interlaced) by entity
2403 * \param[in] tab variable exchanged
2404 * \param[out] local local data
2405 */
2406 /*----------------------------------------------------------------------------*/
2407
2408 void
cs_internal_coupling_exchange_by_cell_id(const cs_internal_coupling_t * cpl,int stride,const cs_real_t tab[],cs_real_t local[])2409 cs_internal_coupling_exchange_by_cell_id(const cs_internal_coupling_t *cpl,
2410 int stride,
2411 const cs_real_t tab[],
2412 cs_real_t local[])
2413 {
2414 int jj;
2415 cs_lnum_t face_id, cell_id;
2416
2417 const cs_lnum_t n_distant = cpl->n_distant;
2418 const cs_lnum_t *faces_distant = cpl->faces_distant;
2419
2420 const cs_mesh_t* m = cs_glob_mesh;
2421
2422 const cs_lnum_t *restrict b_face_cells
2423 = (const cs_lnum_t *restrict)m->b_face_cells;
2424
2425 /* Initialize distant array */
2426
2427 cs_real_t *distant = NULL;
2428 BFT_MALLOC(distant, n_distant*stride, cs_real_t);
2429 for (cs_lnum_t ii = 0; ii < n_distant; ii++) {
2430 face_id = faces_distant[ii];
2431 cell_id = b_face_cells[face_id];
2432 for (jj = 0; jj < stride; jj++)
2433 distant[stride * ii + jj] = tab[stride * cell_id + jj];
2434 }
2435
2436 /* Exchange variable */
2437
2438 cs_internal_coupling_exchange_var(cpl,
2439 stride,
2440 distant,
2441 local);
2442 /* Free memory */
2443 BFT_FREE(distant);
2444 }
2445
2446 /*----------------------------------------------------------------------------*/
2447 /*!
2448 * \brief Exchange variable between groups using face id.
2449 *
2450 * \param[in] cpl pointer to coupling entity
2451 * \param[in] stride number of values (non interlaced) by entity
2452 * \param[in] tab variable exchanged
2453 * \param[out] local local data
2454 */
2455 /*----------------------------------------------------------------------------*/
2456
2457 void
cs_internal_coupling_exchange_by_face_id(const cs_internal_coupling_t * cpl,int stride,const cs_real_t tab[],cs_real_t local[])2458 cs_internal_coupling_exchange_by_face_id(const cs_internal_coupling_t *cpl,
2459 int stride,
2460 const cs_real_t tab[],
2461 cs_real_t local[])
2462 {
2463 const cs_lnum_t n_distant = cpl->n_distant;
2464 const cs_lnum_t *faces_distant = cpl->faces_distant;
2465
2466 /* Initialize distant array */
2467
2468 cs_real_t *distant = NULL;
2469 BFT_MALLOC(distant, n_distant*stride, cs_real_t);
2470 for (cs_lnum_t ii = 0; ii < n_distant; ii++) {
2471 cs_lnum_t face_id = faces_distant[ii];
2472 for (cs_lnum_t jj = 0; jj < stride; jj++)
2473 distant[stride * ii + jj] = tab[stride * face_id + jj];
2474 }
2475
2476 /* Exchange variable */
2477
2478 cs_internal_coupling_exchange_var(cpl,
2479 stride,
2480 distant,
2481 local);
2482 /* Free memory */
2483 BFT_FREE(distant);
2484 }
2485
2486 /*----------------------------------------------------------------------------
2487 * Return pointers to coupling components
2488 *
2489 * parameters:
2490 * cpl <-- pointer to coupling entity
2491 * n_local --> NULL or pointer to component n_local
2492 * faces_local --> NULL or pointer to component faces_local
2493 * n_distant --> NULL or pointer to component n_distant
2494 * faces_distant --> NULL or pointer to component faces_distant
2495 *----------------------------------------------------------------------------*/
2496
2497 void
cs_internal_coupling_coupled_faces(const cs_internal_coupling_t * cpl,cs_lnum_t * n_local,const cs_lnum_t * faces_local[],cs_lnum_t * n_distant,const cs_lnum_t * faces_distant[])2498 cs_internal_coupling_coupled_faces(const cs_internal_coupling_t *cpl,
2499 cs_lnum_t *n_local,
2500 const cs_lnum_t *faces_local[],
2501 cs_lnum_t *n_distant,
2502 const cs_lnum_t *faces_distant[])
2503 {
2504 if (n_local != NULL)
2505 *n_local = cpl->n_local;
2506 if (faces_local != NULL)
2507 *faces_local = cpl->faces_local;
2508 if (n_distant != NULL)
2509 *n_distant = cpl->n_distant;
2510 if (faces_distant != NULL)
2511 *faces_distant = cpl->faces_distant;
2512 }
2513
2514 /*----------------------------------------------------------------------------
2515 * Addition to matrix-vector product in case of internal coupling.
2516 *
2517 * parameters:
2518 * exclude_diag <-- extra diagonal flag
2519 * f <-- associated field pointer
2520 * x <-- vector x in m * x = y
2521 * y <-> vector y in m * x = y
2522 *----------------------------------------------------------------------------*/
2523
2524 void
cs_internal_coupling_spmv_contribution(bool exclude_diag,const cs_field_t * f,const cs_real_t * restrict x,cs_real_t * restrict y)2525 cs_internal_coupling_spmv_contribution(bool exclude_diag,
2526 const cs_field_t *f,
2527 const cs_real_t *restrict x,
2528 cs_real_t *restrict y)
2529 {
2530 cs_lnum_t face_id, cell_id;
2531
2532 const cs_lnum_t *restrict b_face_cells
2533 = (const cs_lnum_t *restrict)cs_glob_mesh->b_face_cells;
2534
2535 int coupling_id = cs_field_get_key_int(f,
2536 cs_field_key_id("coupling_entity"));
2537 const cs_internal_coupling_t *cpl
2538 = cs_internal_coupling_by_id(coupling_id);
2539
2540 const cs_lnum_t n_local = cpl->n_local;
2541 const cs_lnum_t *faces_local = cpl->faces_local;
2542
2543 const int key_cal_opt_id = cs_field_key_id("var_cal_opt");
2544 cs_var_cal_opt_t var_cal_opt;
2545 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
2546 cs_real_t thetap = 0.0;
2547 int idiffp = 0;
2548
2549 if (var_cal_opt.icoupl > 0) {
2550 thetap = var_cal_opt.thetav;
2551 idiffp = var_cal_opt.idiff;
2552 }
2553
2554 /* Exchange x */
2555
2556 cs_real_t *x_j = NULL;
2557 BFT_MALLOC(x_j, f->dim * n_local, cs_real_t);
2558 cs_internal_coupling_exchange_by_cell_id(cpl,
2559 f->dim,
2560 x,
2561 x_j);
2562
2563 /* Compute heq and update y */
2564
2565 cs_real_t *hintp = f->bc_coeffs->hint;
2566 cs_real_t *hextp = f->bc_coeffs->hext;
2567
2568 if (f->dim == 1) {
2569 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
2570 face_id = faces_local[ii];
2571 cell_id = b_face_cells[face_id];
2572
2573 cs_real_t pi = exclude_diag ?
2574 0. : x[cell_id]; /* If exclude_diag, no diagonal term */
2575 cs_real_t pj = x_j[ii];
2576
2577 cs_real_t hint = hintp[face_id];
2578 cs_real_t hext = hextp[face_id];
2579 cs_real_t heq = _calc_heq(hint, hext);
2580
2581 y[cell_id] += thetap * idiffp * heq * (pi - pj);
2582 }
2583
2584 } else if (f->dim == 3) {
2585
2586 cs_real_3_t *_y = (cs_real_3_t *)y;
2587 const cs_real_3_t *_x = (const cs_real_3_t *)x;
2588 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
2589 face_id = faces_local[ii];
2590 cell_id = b_face_cells[face_id];
2591
2592 cs_real_t pi[3];
2593 /* If exclude_diag, no diagonal term */
2594 if (exclude_diag) {
2595 for (cs_lnum_t k = 0; k < 3; k++)
2596 pi[k] = 0.;
2597 } else {
2598 for (cs_lnum_t k = 0; k < 3; k++)
2599 pi[k] = _x[cell_id][k];
2600 }
2601 cs_real_t pj[3] = {x_j[ii], x_j[ii+1], x_j[ii+2]};
2602
2603 cs_real_t hint = hintp[face_id];
2604 cs_real_t hext = hextp[face_id];
2605 cs_real_t heq = _calc_heq(hint, hext);
2606
2607 for (cs_lnum_t k = 0; k < 3; k++)
2608 _y[cell_id][k] += thetap * idiffp * heq * (pi[k] - pj[k]);
2609 }
2610
2611 }
2612
2613 /* Free memory */
2614 BFT_FREE(x_j);
2615 }
2616
2617 /*----------------------------------------------------------------------------
2618 * Add coupling term coordinates to matrix assembler.
2619 *
2620 * parameters:
2621 * coupling_id
2622 * r_g_id <-- global row ids (per cell)
2623 * ma <-> matrix assembler
2624 *----------------------------------------------------------------------------*/
2625
2626 void
cs_internal_coupling_matrix_add_ids(int coupling_id,const cs_gnum_t * r_g_id,cs_matrix_assembler_t * ma)2627 cs_internal_coupling_matrix_add_ids(int coupling_id,
2628 const cs_gnum_t *r_g_id,
2629 cs_matrix_assembler_t *ma)
2630 {
2631 const cs_lnum_t *restrict b_face_cells
2632 = (const cs_lnum_t *restrict)cs_glob_mesh->b_face_cells;
2633 const cs_internal_coupling_t *cpl
2634 = cs_internal_coupling_by_id(coupling_id);
2635
2636 const cs_lnum_t n_distant = cpl->n_distant;
2637 const cs_lnum_t n_local = cpl->n_local;
2638
2639 const cs_lnum_t block_size = 800;
2640 cs_gnum_t g_row_id[800];
2641 cs_gnum_t g_col_id[800];
2642
2643 cs_gnum_t *g_id_l, *g_id_d;
2644 BFT_MALLOC(g_id_l, CS_MAX(n_local, n_distant), cs_gnum_t);
2645 BFT_MALLOC(g_id_d, n_local, cs_gnum_t);
2646
2647 /* local to global preparation and exchange */
2648
2649 for (cs_lnum_t ii = 0; ii < n_distant; ii++) {
2650 cs_lnum_t face_id = cpl->faces_distant[ii];
2651 cs_lnum_t cell_id = b_face_cells[face_id]; /* boundary to cell */
2652 g_id_l[ii] = r_g_id[cell_id];
2653 }
2654
2655 ple_locator_exchange_point_var(cpl->locator,
2656 g_id_l,
2657 g_id_d,
2658 NULL,
2659 sizeof(cs_gnum_t),
2660 1,
2661 0);
2662
2663 /* local side */
2664
2665 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
2666 cs_lnum_t face_id = cpl->faces_local[ii];
2667 cs_lnum_t cell_id = b_face_cells[face_id]; /* boundary to cell */
2668 g_id_l[ii] = r_g_id[cell_id];
2669 }
2670
2671 cs_lnum_t jj = 0;
2672 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
2673 g_row_id[jj] = g_id_l[ii];
2674 g_col_id[jj] = g_id_d[ii];
2675 jj++;
2676 if (jj >= block_size - 1) {
2677 cs_matrix_assembler_add_g_ids(ma, jj, g_row_id, g_col_id);
2678 jj = 0;
2679 }
2680 }
2681 if (jj > 0)
2682 cs_matrix_assembler_add_g_ids(ma, jj, g_row_id, g_col_id);
2683
2684 BFT_FREE(g_id_l);
2685 BFT_FREE(g_id_d);
2686 }
2687
2688 /*----------------------------------------------------------------------------
2689 * Add coupling terms to matrix values assembly.
2690 *
2691 * parameters:
2692 * f <-- associated field
2693 * db_size <-- diagonal block size
2694 * eb_size <-- extra-diagonal block size
2695 * r_g_id <-- global row ids (per cell)
2696 * mav <-> matrix values assembler
2697 *----------------------------------------------------------------------------*/
2698
2699 void
cs_internal_coupling_matrix_add_values(const cs_field_t * f,cs_lnum_t db_size,cs_lnum_t eb_size,const cs_gnum_t r_g_id[],cs_matrix_assembler_values_t * mav)2700 cs_internal_coupling_matrix_add_values(const cs_field_t *f,
2701 cs_lnum_t db_size,
2702 cs_lnum_t eb_size,
2703 const cs_gnum_t r_g_id[],
2704 cs_matrix_assembler_values_t *mav)
2705 {
2706 const cs_lnum_t *restrict b_face_cells
2707 = (const cs_lnum_t *restrict)cs_glob_mesh->b_face_cells;
2708
2709 int coupling_id = cs_field_get_key_int(f, cs_field_key_id("coupling_entity"));
2710 const cs_internal_coupling_t *cpl
2711 = cs_internal_coupling_by_id(coupling_id);
2712
2713 const cs_lnum_t n_distant = cpl->n_distant;
2714 const cs_lnum_t n_local = cpl->n_local;
2715
2716 const int key_cal_opt_id = cs_field_key_id("var_cal_opt");
2717 cs_var_cal_opt_t var_cal_opt;
2718 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
2719 cs_real_t thetap = 0.0;
2720 int idiffp = 0;
2721
2722 if (var_cal_opt.icoupl > 0) {
2723 thetap = var_cal_opt.thetav;
2724 idiffp = var_cal_opt.idiff;
2725 }
2726
2727 /* Compute global ids and exchange coefficient */
2728
2729 cs_real_t *hintp = f->bc_coeffs->hint;
2730 cs_real_t *hextp = f->bc_coeffs->hext;
2731
2732 /* local to global preparation and exchange */
2733
2734 cs_gnum_t *g_id_l, *g_id_d;
2735 BFT_MALLOC(g_id_l, CS_MAX(n_local, n_distant), cs_gnum_t);
2736 BFT_MALLOC(g_id_d, n_local, cs_gnum_t);
2737
2738 for (cs_lnum_t ii = 0; ii < n_distant; ii++) {
2739 cs_lnum_t face_id = cpl->faces_distant[ii];
2740 cs_lnum_t cell_id = b_face_cells[face_id]; /* boundary to cell */
2741 g_id_l[ii] = r_g_id[cell_id];
2742 }
2743
2744 ple_locator_exchange_point_var(cpl->locator,
2745 g_id_l,
2746 g_id_d,
2747 NULL,
2748 sizeof(cs_gnum_t),
2749 1,
2750 0);
2751
2752 /* local side */
2753
2754 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
2755 cs_lnum_t face_id = cpl->faces_local[ii];
2756 cs_lnum_t cell_id = b_face_cells[face_id]; /* boundary to cell */
2757 g_id_l[ii] = r_g_id[cell_id];
2758 }
2759
2760 /* Compute contributions to matrix */
2761
2762 const cs_lnum_t block_size = 514;
2763 cs_gnum_t d_g_row_id[514];
2764 cs_real_t d_aij[514];
2765 cs_gnum_t e_g_row_id[514];
2766 cs_gnum_t e_g_col_id[514];
2767 cs_real_t e_aij[514];
2768
2769 cs_lnum_t db_stride = db_size*db_size;
2770 cs_lnum_t eb_stride = eb_size*eb_size;
2771
2772 assert(block_size > db_stride && block_size > eb_size);
2773
2774 cs_lnum_t jj = 0, kk = 0, db_fill = 0, eb_fill = 0;
2775 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
2776 cs_lnum_t face_id = cpl->faces_local[ii];
2777
2778 cs_real_t hint = hintp[face_id];
2779 cs_real_t hext = hextp[face_id];
2780 cs_real_t heq = _calc_heq(hint, hext);
2781 cs_real_t c = thetap * idiffp * heq;
2782
2783 d_g_row_id[jj] = g_id_l[ii];
2784 e_g_row_id[kk] = g_id_l[ii]; e_g_col_id[kk] = g_id_d[ii];
2785
2786 for (cs_lnum_t ib = 0; ib < db_stride; ib++)
2787 d_aij[db_fill + ib] = 0;
2788 for (cs_lnum_t ib = 0; ib < db_size; ib++)
2789 d_aij[db_fill + ib*(db_size + 1)] = c;
2790
2791 for (cs_lnum_t ib = 0; ib < eb_stride; ib++)
2792 e_aij[eb_fill + ib] = 0;
2793 for (cs_lnum_t ib = 0; ib < eb_size; ib++)
2794 e_aij[eb_fill + ib*(eb_size + 1)] = -c;
2795
2796 jj += 1;
2797 kk += 1;
2798 db_fill += db_stride;
2799 eb_fill += eb_stride;
2800
2801 if (db_fill >= block_size - 1) {
2802 cs_matrix_assembler_values_add_g(mav, jj, d_g_row_id, d_g_row_id, d_aij);
2803 jj = 0;
2804 db_fill = 0;
2805 }
2806
2807 if (eb_fill >= block_size - 1) {
2808 cs_matrix_assembler_values_add_g(mav, kk, e_g_row_id, e_g_col_id, e_aij);
2809 kk = 0;
2810 eb_fill = 0;
2811 }
2812 }
2813
2814 cs_matrix_assembler_values_add_g(mav, jj, d_g_row_id, d_g_row_id, d_aij);
2815 cs_matrix_assembler_values_add_g(mav, kk, e_g_row_id, e_g_col_id, e_aij);
2816
2817 /* Free memory */
2818
2819 BFT_FREE(g_id_l);
2820 BFT_FREE(g_id_d);
2821 }
2822
2823 /*----------------------------------------------------------------------------*/
2824 /*!
2825 * \brief Setup internal coupling related parameters.
2826 */
2827 /*----------------------------------------------------------------------------*/
2828
2829 void
cs_internal_coupling_setup(void)2830 cs_internal_coupling_setup(void)
2831 {
2832 /* Call deprecated functions first */
2833 cs_user_internal_coupling_add_volumes(cs_glob_mesh);
2834 cs_user_internal_coupling_from_disjoint_meshes(cs_glob_mesh);
2835
2836 /* Now do setup proper */
2837
2838 if (_n_internal_couplings < 1)
2839 return;
2840
2841 int field_id;
2842 cs_field_t *f;
2843 const int coupling_key_id = cs_field_key_id("coupling_entity");
2844 int coupling_id = 0;
2845
2846 const int key_cal_opt_id = cs_field_key_id("var_cal_opt");
2847 cs_var_cal_opt_t var_cal_opt;
2848
2849 const int n_fields = cs_field_n_fields();
2850
2851 /* Definition of coupling_ids as keys of variable fields */
2852
2853 coupling_id = 0;
2854 for (field_id = 0; field_id < n_fields; field_id++) {
2855 f = cs_field_by_id(field_id);
2856 if (f->type & CS_FIELD_VARIABLE) {
2857 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
2858 if (var_cal_opt.icoupl > 0) {
2859 cs_field_set_key_int(f, coupling_key_id, coupling_id);
2860 // coupling_id++;
2861 }
2862 }
2863 }
2864
2865 /* Initialization of coupling entities */
2866
2867 coupling_id = 0;
2868 for (field_id = 0; field_id < n_fields; field_id++) {
2869 f = cs_field_by_id(field_id);
2870 if (f->type & CS_FIELD_VARIABLE) {
2871 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
2872 if (var_cal_opt.icoupl > 0) {
2873 coupling_id++;
2874 }
2875 }
2876 }
2877 }
2878
2879 /*----------------------------------------------------------------------------*/
2880 /*!
2881 * \brief Initialize internal coupling related structures.
2882 */
2883 /*----------------------------------------------------------------------------*/
2884
2885 void
cs_internal_coupling_initialize(void)2886 cs_internal_coupling_initialize(void)
2887 {
2888 for (int i = 0; i < _n_internal_couplings; i++) {
2889 cs_internal_coupling_t *cpl = _internal_coupling + i;
2890 _locator_initialize(cs_glob_mesh, cpl);
2891 _initialize_coupled_faces(cpl);
2892 }
2893 }
2894
2895 /*----------------------------------------------------------------------------
2896 * Log information about a given internal coupling entity
2897 *
2898 * parameters:
2899 * cpl <-- pointer to coupling entity
2900 *----------------------------------------------------------------------------*/
2901
2902 void
cs_internal_coupling_log(const cs_internal_coupling_t * cpl)2903 cs_internal_coupling_log(const cs_internal_coupling_t *cpl)
2904 {
2905 if (cpl == NULL)
2906 return;
2907
2908 cs_gnum_t n_local = cpl->n_local;
2909
2910 cs_parall_counter(&n_local, 1);
2911
2912 if (cpl->cells_criteria != NULL)
2913 bft_printf(" Cell group selection criterion: %s\n",
2914 cpl->cells_criteria);
2915
2916 if (cpl->faces_criteria != NULL)
2917 bft_printf(" Face group selection criterion: %s\n",
2918 cpl->faces_criteria);
2919
2920 if (cpl->interior_faces_group_name != NULL)
2921 bft_printf(" Assign interior faces group name: %s\n",
2922 cpl->interior_faces_group_name);
2923
2924 if (cpl->exterior_faces_group_name != NULL)
2925 bft_printf(" Assign interior faces group name: %s\n",
2926 cpl->exterior_faces_group_name);
2927
2928 if (cpl->n_volume_zones > 0) {
2929 bft_printf(" Volume zones:\n");
2930 for (int i = 0; i < cpl->n_volume_zones; i++) {
2931 const cs_zone_t *z = cs_volume_zone_by_id(cpl->volume_zone_ids[i]);
2932 bft_printf(" %s\n", z->name);
2933 }
2934 }
2935
2936 bft_printf("\n"
2937 " Locator: n dist points (total coupled boundary faces) = %llu\n",
2938 (unsigned long long)n_local);
2939 }
2940
2941 /*----------------------------------------------------------------------------
2942 * Print informations about all coupling entities
2943 *
2944 * parameters:
2945 * cpl <-- pointer to coupling entity
2946 *----------------------------------------------------------------------------*/
2947
2948 void
cs_internal_coupling_dump(void)2949 cs_internal_coupling_dump(void)
2950 {
2951 cs_internal_coupling_t* cpl;
2952
2953 if (_n_internal_couplings == 0)
2954 return;
2955
2956 bft_printf("\n Internal coupling\n");
2957 for (int cpl_id = 0; cpl_id < _n_internal_couplings; cpl_id++) {
2958 cpl = _internal_coupling + cpl_id;
2959 bft_printf(" coupling_id = %d\n", cpl_id);
2960 cs_internal_coupling_log(cpl);
2961 }
2962 }
2963
2964 /*----------------------------------------------------------------------------
2965 * Add preprocessing operations required by coupling volume using given
2966 * criteria.
2967 *
2968 * The volume is separated from the rest of the domain with inserted
2969 * boundaries.
2970 *
2971 * parameters:
2972 * mesh <-> pointer to mesh structure to modify
2973 *----------------------------------------------------------------------------*/
2974
2975 void
cs_internal_coupling_preprocess(cs_mesh_t * mesh)2976 cs_internal_coupling_preprocess(cs_mesh_t *mesh)
2977 {
2978 for (int i = 0; i < _n_internal_couplings; i++) {
2979 cs_internal_coupling_t *cpl = _internal_coupling + i;
2980 if ( (cpl->cells_criteria != NULL || cpl->n_volume_zones > 0)
2981 && cpl->faces_criteria == NULL) {
2982 /* Insert wall boundaries,
2983 * locators are initialized afterwards */
2984 _volume_initialize_insert_boundary(mesh, cpl);
2985 }
2986 }
2987 }
2988
2989 /*----------------------------------------------------------------------------
2990 * Define face to face mappings for internal couplings.
2991 *
2992 * parameters:
2993 * mesh <-> pointer to mesh structure to modify
2994 *----------------------------------------------------------------------------*/
2995
2996 void
cs_internal_coupling_map(cs_mesh_t * mesh)2997 cs_internal_coupling_map(cs_mesh_t *mesh)
2998 {
2999 /* Initialization of locators for all coupling entities */
3000
3001 for (int cpl_id = 0; cpl_id < _n_internal_couplings; cpl_id++) {
3002 cs_internal_coupling_t *cpl = _internal_coupling + cpl_id;
3003 if (cpl->faces_criteria == NULL)
3004 _auto_group_name(cpl, cpl_id);
3005 _volume_face_initialize(mesh, cpl);
3006 }
3007 }
3008
3009 /*----------------------------------------------------------------------------
3010 * Define coupling entity using given criteria.
3011 *
3012 * parameters:
3013 * f_id <-- id of the field
3014 *----------------------------------------------------------------------------*/
3015
3016 void
cs_internal_coupling_add_entity(int f_id)3017 cs_internal_coupling_add_entity(int f_id)
3018 {
3019 const int key_cal_opt_id = cs_field_key_id("var_cal_opt");
3020 cs_var_cal_opt_t var_cal_opt;
3021
3022 cs_field_t* f = cs_field_by_id(f_id);
3023
3024 if (f->type & CS_FIELD_VARIABLE) {
3025 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
3026 var_cal_opt.icoupl = 1;
3027 cs_field_set_key_struct(f, key_cal_opt_id, &var_cal_opt);
3028 }
3029 else
3030 bft_error(__FILE__, __LINE__, 0,
3031 "field id = %d provided is invalid."
3032 " The field must be a variable.",
3033 f_id);
3034 }
3035
3036 /*----------------------------------------------------------------------------*/
3037 /*!
3038 * \brief Update internal coupling coefficients of the field of the
3039 * given id using given boundary exchange coefficients passed by face id.
3040 *
3041 * \param[in] field_id field id
3042 * \param[in] hbnd boundary exchange coefficients passed by face id
3043 */
3044 /*----------------------------------------------------------------------------*/
3045
3046 void
cs_ic_field_set_exchcoeff(const int field_id,const cs_real_t * hbnd)3047 cs_ic_field_set_exchcoeff(const int field_id,
3048 const cs_real_t *hbnd)
3049 {
3050 const cs_real_t* b_face_surf = cs_glob_mesh_quantities->b_face_surf;
3051
3052 const cs_field_t* f = cs_field_by_id(field_id);
3053
3054 const int coupling_key_id = cs_field_key_id("coupling_entity");
3055 int coupling_id = cs_field_get_key_int(f,
3056 coupling_key_id);
3057 const cs_internal_coupling_t *cpl
3058 = cs_internal_coupling_by_id(coupling_id);
3059
3060 const cs_lnum_t n_local = cpl->n_local;
3061 const cs_lnum_t *faces_local = cpl->faces_local;
3062 cs_real_t *hint = f->bc_coeffs->hint;
3063 cs_real_t *hext = f->bc_coeffs->hext;
3064
3065 cs_real_t *hextloc = NULL;
3066 BFT_MALLOC(hextloc, n_local, cs_real_t);
3067
3068 /* Exchange hbnd */
3069 cs_internal_coupling_exchange_by_face_id(cpl,
3070 1, /* Dimension */
3071 hbnd,
3072 hextloc);
3073
3074 /* Compute hint and hext */
3075 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
3076 cs_lnum_t face_id = faces_local[ii];
3077 cs_real_t surf = b_face_surf[face_id];
3078 hint[face_id] = hbnd[face_id] * surf;
3079 hext[face_id] = hextloc[ii] * surf;
3080 }
3081
3082 BFT_FREE(hextloc);
3083 }
3084
3085 /*----------------------------------------------------------------------------*/
3086 /*!
3087 * \brief Get distant data using face id at all coupling faces for a given
3088 * field id.
3089 *
3090 * \param[in] field_id field id
3091 * \param[in] stride number of values (interlaced) by entity
3092 * \param[in] tab_distant exchanged data by face id
3093 * \param[out] tab_local local data by face id
3094 */
3095 /*----------------------------------------------------------------------------*/
3096
3097 void
cs_ic_field_dist_data_by_face_id(const int field_id,int stride,const cs_real_t tab_distant[],cs_real_t tab_local[])3098 cs_ic_field_dist_data_by_face_id(const int field_id,
3099 int stride,
3100 const cs_real_t tab_distant[],
3101 cs_real_t tab_local[])
3102 {
3103 const cs_field_t* f = cs_field_by_id(field_id);
3104
3105 const int coupling_key_id = cs_field_key_id("coupling_entity");
3106 int coupling_id = cs_field_get_key_int(f,
3107 coupling_key_id);
3108 const cs_internal_coupling_t *cpl
3109 = cs_internal_coupling_by_id(coupling_id);
3110 const cs_lnum_t n_local = cpl->n_local;
3111 const cs_lnum_t *faces_local = cpl->faces_local;
3112
3113 cs_real_t *local = NULL;
3114 BFT_MALLOC(local, n_local, cs_real_t);
3115
3116 /* Exchange distant data by face id */
3117 cs_internal_coupling_exchange_by_face_id(cpl,
3118 stride,
3119 tab_distant,
3120 local);
3121
3122 /* Save by face id */
3123 for (cs_lnum_t ii = 0; ii < n_local; ii++) {
3124 cs_lnum_t face_id = faces_local[ii];
3125 for (cs_lnum_t jj = 0; jj < stride; jj++)
3126 tab_local[stride * face_id + jj] = local[stride * ii + jj];
3127 }
3128
3129 BFT_FREE(local);
3130 }
3131
3132 /*----------------------------------------------------------------------------*/
3133 /*!
3134 * \brief Define volumes as internal coupling zones.
3135 *
3136 * These zones will be separated from the rest of the domain using automatically
3137 * defined thin walls.
3138 *
3139 * \deprecated
3140 * move contents to\ref cs_user_internal_coupling instead.
3141 *
3142 * \param[in, out] mesh pointer to a cs_mesh_t structure
3143 */
3144 /*----------------------------------------------------------------------------*/
3145
3146 void
cs_user_internal_coupling_add_volumes(cs_mesh_t * mesh)3147 cs_user_internal_coupling_add_volumes(cs_mesh_t *mesh)
3148 {
3149 CS_UNUSED(mesh);
3150 }
3151
3152 /*----------------------------------------------------------------------------*/
3153 /*!
3154 * \brief Define volumes from separated meshes as internal coupling zones.
3155 *
3156 * These zones must be disjoint and the face selection criteria must be specified.
3157 *
3158 * \deprecated
3159 * move contents to\ref cs_user_internal_coupling instead.
3160 *
3161 * \param[in, out] mesh pointer to a cs_mesh_t structure
3162 */
3163 /*----------------------------------------------------------------------------*/
3164
3165 void
cs_user_internal_coupling_from_disjoint_meshes(cs_mesh_t * mesh)3166 cs_user_internal_coupling_from_disjoint_meshes(cs_mesh_t *mesh)
3167 {
3168 CS_UNUSED(mesh);
3169 }
3170
3171 /*----------------------------------------------------------------------------*/
3172
3173 END_C_DECLS
3174