1 /*============================================================================
2 * Boundary condition handling.
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <stdarg.h>
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <string.h>
37 #include <assert.h>
38 #include <math.h>
39
40 #if defined(HAVE_MPI)
41 #include <mpi.h>
42 #endif
43
44 /*----------------------------------------------------------------------------
45 * Local headers
46 *----------------------------------------------------------------------------*/
47
48 #include <ple_locator.h>
49
50 #include "bft_mem.h"
51 #include "bft_error.h"
52 #include "bft_printf.h"
53
54 #include "cs_base.h"
55 #include "cs_boundary.h"
56 #include "cs_coupling.h"
57 #include "cs_gradient.h"
58 #include "cs_gui_util.h"
59 #include "cs_field.h"
60 #include "cs_field_default.h"
61 #include "cs_field_operator.h"
62 #include "cs_field_pointer.h"
63 #include "cs_flag_check.h"
64 #include "cs_halo.h"
65 #include "cs_math.h"
66 #include "cs_mesh.h"
67 #include "cs_mesh_connect.h"
68 #include "cs_mesh_quantities.h"
69 #include "cs_parall.h"
70 #include "cs_parameters.h"
71 #include "cs_physical_model.h"
72 #include "cs_prototypes.h"
73 #include "cs_post.h"
74 #include "cs_xdef_eval_at_zone.h"
75 #include "fvm_nodal.h"
76
77 /*----------------------------------------------------------------------------
78 * Header for the current file
79 *----------------------------------------------------------------------------*/
80
81 #include "cs_boundary_conditions.h"
82
83 /*----------------------------------------------------------------------------*/
84
85 BEGIN_C_DECLS
86
87 /*=============================================================================
88 * Additional doxygen documentation
89 *============================================================================*/
90
91 /*!
92 \file cs_boundary_conditions.c
93 Boundary condition handling.
94 */
95
96 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
97
98 /*=============================================================================
99 * Local Macro Definitions
100 *============================================================================*/
101
102 /*=============================================================================
103 * Local Type Definitions
104 *============================================================================*/
105
106 typedef struct {
107
108 int bc_location_id; /* location id of boundary zone */
109 int source_location_id; /* location id of source elements */
110 cs_real_t coord_shift[3]; /* coordinates shift relative to
111 selected boundary faces */
112 double tolerance; /* search tolerance */
113
114 ple_locator_t *locator; /* associated locator */
115
116 } cs_bc_map_t;
117
118 /*----------------------------------------------------------------------------
119 * Local Structure Definitions
120 *----------------------------------------------------------------------------*/
121
122 /*============================================================================
123 * Global variables
124 *============================================================================*/
125
126 static int *_bc_type;
127
128 const int *cs_glob_bc_type = NULL;
129
130 static int *_bc_face_zone;
131
132 const int *cs_glob_bc_face_zone = NULL;
133
134 static int _n_bc_maps = 0;
135 static cs_bc_map_t *_bc_maps = NULL;
136
137 /*============================================================================
138 * Prototypes for functions intended for use only by Fortran wrappers.
139 * (descriptions follow, with function bodies).
140 *============================================================================*/
141
142 void
143 cs_f_boundary_conditions_mapped_set(int field_id,
144 ple_locator_t *locator,
145 cs_mesh_location_type_t location_type,
146 int enforce_balance,
147 int interpolate,
148 cs_lnum_t n_faces,
149 const cs_lnum_t *faces,
150 cs_real_t *balance_w,
151 int nvar,
152 cs_real_t *rcodcl);
153
154 void
155 cs_f_boundary_conditions_get_pointers(int **itypfb,
156 int **izfppp);
157
158 /*============================================================================
159 * Private function definitions
160 *============================================================================*/
161
162 /*----------------------------------------------------------------------------*/
163 /*!
164 * \brief Build or update map of shifted boundary face coordinates on
165 * cells or boundary faces for automatic interpolation.
166 *
167 * \warning
168 *
169 * This function does not currently update the location in case of moving
170 * meshes. It could be modifed to do so.
171 *
172 * \param[in] map_id id of defined map
173 */
174 /*----------------------------------------------------------------------------*/
175
176 static void
_update_bc_map(int map_id)177 _update_bc_map(int map_id)
178 {
179 assert(map_id > -1 && map_id < _n_bc_maps);
180
181 cs_bc_map_t *bc_map = _bc_maps + map_id;
182
183 if (bc_map->locator != NULL)
184 return;
185
186 cs_mesh_location_type_t location_type
187 = cs_mesh_location_get_type(bc_map->source_location_id);
188
189 cs_lnum_t n_location_elts
190 = cs_mesh_location_get_n_elts(bc_map->source_location_id)[0];
191 cs_lnum_t n_faces
192 = cs_mesh_location_get_n_elts(bc_map->bc_location_id)[0];
193
194 const cs_lnum_t *location_elts
195 = cs_mesh_location_get_elt_ids_try(bc_map->source_location_id);
196 const cs_lnum_t *faces
197 = cs_mesh_location_get_elt_ids_try(bc_map->bc_location_id);
198
199 bc_map->locator = cs_boundary_conditions_map(location_type,
200 n_location_elts,
201 n_faces,
202 location_elts,
203 faces,
204 &(bc_map->coord_shift),
205 0,
206 bc_map->tolerance);
207 }
208
209 /*----------------------------------------------------------------------------
210 * Compute balance at inlet
211 *
212 * parameters:
213 * var_id <-- variable id
214 * f <-- associated field
215 * m <-- associated mesh
216 * mq <-- mesh quantities
217 * enforce_balance <-- balance handling option:
218 * 0: values are simply mapped
219 * 1: values are mapped, then multiplied
220 * by a constant factor so that their
221 * surface integral on selected faces
222 * is preserved (relative to the
223 * input values)
224 * 2: as 1, but with a boundary-defined
225 * weight, defined by balance_w
226 * 3: as 1, but with a cell-defined
227 * weight, defined by balance_w
228 * n_faces <-- number of selected boundary faces
229 * faces <-- list of selected boundary faces (0 to n-1),
230 * or NULL if no indirection is needed
231 * balance_w <-- optional balance weight, or NULL
232 * nvar <-- number of variables requiring BC's
233 * rcodcl <-- boundary condition values
234 * inlet_sum --> inlet sum
235 *----------------------------------------------------------------------------*/
236
237 static void
_inlet_sum(int var_id,const cs_field_t * f,const cs_mesh_t * m,const cs_mesh_quantities_t * mq,int enforce_balance,cs_lnum_t n_faces,const cs_lnum_t * faces,cs_real_t * balance_w,cs_lnum_t nvar,cs_real_t rcodcl[],cs_real_t inlet_sum[])238 _inlet_sum(int var_id,
239 const cs_field_t *f,
240 const cs_mesh_t *m,
241 const cs_mesh_quantities_t *mq,
242 int enforce_balance,
243 cs_lnum_t n_faces,
244 const cs_lnum_t *faces,
245 cs_real_t *balance_w,
246 cs_lnum_t nvar,
247 cs_real_t rcodcl[],
248 cs_real_t inlet_sum[])
249 {
250 CS_UNUSED(nvar);
251
252 const int dim = f->dim;
253 const cs_lnum_t n_b_faces = m->n_b_faces;
254 const cs_real_t *f_surf = mq->b_f_face_surf;
255
256 /* Get field's variable id */
257
258 assert(dim <= 9);
259
260 for (cs_lnum_t j = 0; j < dim; j++) {
261
262 cs_real_t *_rcodcl = rcodcl + (var_id+j)*n_b_faces;
263
264 inlet_sum[j] = 0.;
265
266 if (enforce_balance == 1) {
267 for (cs_lnum_t i = 0; i < n_faces; i++) {
268 const cs_lnum_t f_id = (faces != NULL) ? faces[i] : i;
269 inlet_sum[j] +=_rcodcl[f_id]*f_surf[f_id];
270 }
271 }
272 else if (enforce_balance == 2) {
273 for (cs_lnum_t i = 0; i < n_faces; i++) {
274 const cs_lnum_t f_id = (faces != NULL) ? faces[i] : i;
275 inlet_sum[j] +=_rcodcl[f_id]*f_surf[f_id]*balance_w[f_id];
276 }
277 }
278 else if (enforce_balance == 3) {
279 for (cs_lnum_t i = 0; i < n_faces; i++) {
280 const cs_lnum_t f_id = (faces != NULL) ? faces[i] : i;
281 const cs_lnum_t c_id = m->b_face_cells[f_id];
282 inlet_sum[j] +=_rcodcl[f_id]*f_surf[f_id]*balance_w[c_id];
283 }
284 }
285
286 }
287
288 cs_parall_sum(dim, CS_REAL_TYPE, inlet_sum);
289 }
290
291 /*----------------------------------------------------------------------------*/
292 /*!
293 * \brief Compute the values of the homogeneous Dirichlet BCs
294 *
295 * \param[in] mesh pointer to mesh structure
296 * \param[in] boundaries pointer to associated boundaries
297 * \param[in] eqp pointer to a cs_equation_param_t
298 * \param[in] def pointer to a boundary condition definition
299 * \param[in] var_id matching variable id
300 * \param[in] icodcl_m multiplier for assigned icodcl values (1 or -1)
301 * \param[in, out] icodcl boundary conditions type array
302 * \param[in, out] rcodcl boundary conditions values array
303 */
304 /*----------------------------------------------------------------------------*/
305
306 static void
_compute_hmg_dirichlet_bc(const cs_mesh_t * mesh,const cs_boundary_t * boundaries,const cs_equation_param_t * eqp,const cs_xdef_t * def,cs_lnum_t var_id,int icodcl_m,int * icodcl,double * rcodcl)307 _compute_hmg_dirichlet_bc(const cs_mesh_t *mesh,
308 const cs_boundary_t *boundaries,
309 const cs_equation_param_t *eqp,
310 const cs_xdef_t *def,
311 cs_lnum_t var_id,
312 int icodcl_m,
313 int *icodcl,
314 double *rcodcl)
315 {
316 assert(eqp->dim == def->dim);
317
318 const cs_lnum_t n_b_faces = mesh->n_b_faces;
319 const cs_zone_t *bz = cs_boundary_zone_by_id(def->z_id);
320 const cs_lnum_t *face_ids = bz->elt_ids;
321
322 cs_boundary_type_t boundary_type
323 = boundaries->types[cs_boundary_id_by_zone_id(boundaries, def->z_id)];
324
325 int bc_type = (boundary_type & CS_BOUNDARY_WALL) ? 5 : 1;
326 if (boundary_type & CS_BOUNDARY_ROUGH_WALL)
327 bc_type = 6;
328
329 bc_type *= icodcl_m;
330
331 for (int coo_id = 0; coo_id < def->dim; coo_id++) {
332
333 int *_icodcl = icodcl + (var_id+coo_id)*n_b_faces;
334 cs_real_t *_rcodcl1 = rcodcl + (var_id+coo_id)*n_b_faces;
335
336 # pragma omp parallel for if (bz->n_elts > CS_THR_MIN)
337 for (cs_lnum_t i = 0; i < bz->n_elts; i++) {
338 const cs_lnum_t elt_id = (face_ids == NULL) ? i : face_ids[i];
339 _icodcl[elt_id] = bc_type;
340 _rcodcl1[elt_id] = 0;
341 }
342 }
343 }
344
345 /*----------------------------------------------------------------------------*/
346 /*!
347 * \brief Compute the values of the Dirichlet BCs
348 *
349 * \param[in] mesh pointer to mesh structure
350 * \param[in] boundaries pointer to associated boundaries
351 * \param[in] f pointer to associated field
352 * \param[in] eqp pointer to a cs_equation_param_t
353 * \param[in] def pointer to a boundary condition definition
354 * \param[in] description description string (for error logging)
355 * \param[in] var_id matching variable id
356 * \param[in] t_eval time at which one evaluates the boundary cond.
357 * \param[in] icodcl_m multiplier for assigned icodcl values (1 or -1)
358 * \param[in, out] icodcl boundary conditions type array
359 * \param[in, out] rcodcl boundary conditions values array
360 * \param eval_buf evaluation bufferb (work array)
361 */
362 /*----------------------------------------------------------------------------*/
363
364 static void
_compute_dirichlet_bc(const cs_mesh_t * mesh,const cs_boundary_t * boundaries,const cs_field_t * f,const cs_equation_param_t * eqp,const cs_xdef_t * def,const char * description,cs_lnum_t var_id,cs_real_t t_eval,int icodcl_m,int * icodcl,double * rcodcl,cs_real_t eval_buf[])365 _compute_dirichlet_bc(const cs_mesh_t *mesh,
366 const cs_boundary_t *boundaries,
367 const cs_field_t *f,
368 const cs_equation_param_t *eqp,
369 const cs_xdef_t *def,
370 const char *description,
371 cs_lnum_t var_id,
372 cs_real_t t_eval,
373 int icodcl_m,
374 int *icodcl,
375 double *rcodcl,
376 cs_real_t eval_buf[])
377 {
378 CS_UNUSED(t_eval);
379
380 const cs_lnum_t n_b_faces = mesh->n_b_faces;
381 const cs_zone_t *bz = cs_boundary_zone_by_id(def->z_id);
382 const cs_lnum_t *elt_ids = bz->elt_ids;
383 const cs_lnum_t n_elts = bz->n_elts;
384 const cs_lnum_t def_dim = def->dim;
385
386 assert(eqp->dim == def->dim);
387
388 cs_boundary_type_t boundary_type
389 = boundaries->types[cs_boundary_id_by_zone_id(boundaries, def->z_id)];
390
391 int bc_type = (boundary_type & CS_BOUNDARY_WALL) ? 5 : 1;
392 if (boundary_type & CS_BOUNDARY_ROUGH_WALL)
393 bc_type = 6;
394
395 bc_type *= icodcl_m;
396
397 if (f->dim != def_dim) {
398 bft_error(__FILE__, __LINE__, 0,
399 _(" %s: Boundary condition definition:\n"
400 " field %s, zone %s, type %s;\n"
401 " dimension %d does not match field dimension (%d)."),
402 __func__, f->name, bz->name, cs_xdef_type_get_name(def->type),
403 def_dim, f->dim);
404 }
405
406 switch(def->type) {
407
408 case CS_XDEF_BY_VALUE: /* direct application (optimization)
409 for constant value */
410 {
411 const cs_real_t *constant_val = (cs_real_t *)def->context;
412
413 for (cs_lnum_t coo_id = 0; coo_id < def_dim; coo_id++) {
414
415 int *_icodcl = icodcl + (var_id+coo_id)*n_b_faces;
416 cs_real_t *_rcodcl1 = rcodcl + (var_id+coo_id)*n_b_faces;
417
418 # pragma omp parallel for if (n_elts > CS_THR_MIN)
419 for (cs_lnum_t i = 0; i < n_elts; i++) {
420 const cs_lnum_t elt_id = elt_ids[i];
421 _icodcl[elt_id] = bc_type;
422 _rcodcl1[elt_id] = constant_val[coo_id];
423 }
424
425 }
426 }
427 break;
428
429 default:
430 {
431 cs_xdef_eval_at_zone(def,
432 description,
433 t_eval,
434 true, /* dense */
435 eval_buf);
436
437 const cs_lnum_t _dim = def->dim;
438
439 for (cs_lnum_t coo_id = 0; coo_id < _dim; coo_id++) {
440
441 int *_icodcl = icodcl + (var_id+coo_id)*n_b_faces;
442 cs_real_t *_rcodcl1 = rcodcl + (var_id+coo_id)*n_b_faces;
443
444 # pragma omp parallel for if (n_elts > CS_THR_MIN)
445 for (cs_lnum_t i = 0; i < n_elts; i++) {
446 const cs_lnum_t elt_id = elt_ids[i];
447 _icodcl[elt_id] = bc_type;
448 _rcodcl1[elt_id] = eval_buf[_dim*i + coo_id];
449 }
450
451 }
452 }
453 break;
454 }
455 }
456
457 /*----------------------------------------------------------------------------*/
458 /*!
459 * \brief Compute the values of the homogeneous Neumann BCs
460 *
461 * \param[in] mesh pointer to mesh structure
462 * \param[in] eqp pointer to a cs_equation_param_t
463 * \param[in] def pointer to a boundary condition definition
464 * \param[in] nvar number of variables
465 * \param[in] var_id matching variable id
466 * \param[in, out] icodcl boundary conditions type array
467 * \param[in, out] rcodcl boundary conditions values array
468 */
469 /*----------------------------------------------------------------------------*/
470
471 static void
_compute_hmg_neumann_bc(const cs_mesh_t * mesh,const cs_xdef_t * def,cs_lnum_t nvar,cs_lnum_t var_id,int * icodcl,double * rcodcl)472 _compute_hmg_neumann_bc(const cs_mesh_t *mesh,
473 const cs_xdef_t *def,
474 cs_lnum_t nvar,
475 cs_lnum_t var_id,
476 int *icodcl,
477 double *rcodcl)
478 {
479 const cs_lnum_t n_b_faces = mesh->n_b_faces;
480 const cs_zone_t *bz = cs_boundary_zone_by_id(def->z_id);
481 const cs_lnum_t *elt_ids = bz->elt_ids;
482 const cs_lnum_t n_elts = bz->n_elts;
483
484 for (int coo_id = 0; coo_id < def->dim; coo_id++) {
485
486 int *_icodcl = icodcl + (var_id+coo_id)*n_b_faces;
487 cs_real_t *_rcodcl3 = rcodcl + 2*n_b_faces*nvar
488 + (var_id+coo_id)*n_b_faces;
489
490 # pragma omp parallel for if (n_elts > CS_THR_MIN)
491 for (cs_lnum_t i = 0; i < n_elts; i++) {
492 const cs_lnum_t elt_id = elt_ids[i];
493 _icodcl[elt_id] = 3;
494 _rcodcl3[elt_id] = 0;
495 }
496 }
497 }
498
499 /*----------------------------------------------------------------------------*/
500 /*!
501 * \brief Compute the values of the Neumann BCs
502 *
503 * \param[in] mesh pointer to mesh structure
504 * \param[in] eqp pointer to a cs_equation_param_t
505 * \param[in] def pointer to a boundary condition definition
506 * \param[in] description description string (for error logging)
507 * \param[in] nvar number of variables
508 * \param[in] var_id matching variable id
509 * \param[in] t_eval time at which one evaluates the boundary cond.
510 * \param[in, out] icodcl boundary conditions type array
511 * \param[in, out] rcodcl boundary conditions values array
512 * \param eval_buf evaluation bufferb (work array)
513 */
514 /*----------------------------------------------------------------------------*/
515
516 static void
_compute_neumann_bc(const cs_mesh_t * mesh,const cs_equation_param_t * eqp,const cs_xdef_t * def,const char * description,cs_lnum_t nvar,cs_lnum_t var_id,cs_real_t t_eval,int * icodcl,double * rcodcl,cs_real_t eval_buf[])517 _compute_neumann_bc(const cs_mesh_t *mesh,
518 const cs_equation_param_t *eqp,
519 const cs_xdef_t *def,
520 const char *description,
521 cs_lnum_t nvar,
522 cs_lnum_t var_id,
523 cs_real_t t_eval,
524 int *icodcl,
525 double *rcodcl,
526 cs_real_t eval_buf[])
527 {
528 const cs_lnum_t n_b_faces = mesh->n_b_faces;
529 const cs_zone_t *bz = cs_boundary_zone_by_id(def->z_id);
530 const cs_lnum_t *elt_ids = bz->elt_ids;
531 const cs_lnum_t n_elts = bz->n_elts;
532 const cs_lnum_t def_dim = def->dim;
533
534 assert(eqp->dim == def->dim);
535
536 const int bc_type = 3;
537
538 switch(def->type) {
539
540 case CS_XDEF_BY_VALUE:
541 {
542 const int dim = def->dim;
543 const cs_real_t *constant_vals = (cs_real_t *)def->context;
544
545 for (int coo_id = 0; coo_id < dim; coo_id++) {
546
547 const cs_real_t value = constant_vals[coo_id];
548
549 int *_icodcl = icodcl + (var_id+coo_id)*n_b_faces;
550 cs_real_t *_rcodcl3 = rcodcl + 2*n_b_faces*nvar
551 + (var_id+coo_id)*n_b_faces;
552
553 # pragma omp parallel for if (bz->n_elts > CS_THR_MIN)
554 for (cs_lnum_t i = 0; i < bz->n_elts; i++) {
555 const cs_lnum_t elt_id = elt_ids[i];
556 _icodcl[elt_id] = bc_type;
557 _rcodcl3[elt_id] = value;
558 }
559
560 }
561 }
562 break;
563
564 default:
565 {
566 cs_xdef_eval_at_zone(def,
567 description,
568 t_eval,
569 true, /* dense */
570 eval_buf);
571
572 const cs_lnum_t _dim = def->dim;
573
574 for (cs_lnum_t coo_id = 0; coo_id < _dim; coo_id++) {
575
576 int *_icodcl = icodcl + (var_id+coo_id)*n_b_faces;
577 cs_real_t *_rcodcl3 = rcodcl + 2*n_b_faces*nvar
578 + (var_id+coo_id)*n_b_faces;
579
580 # pragma omp parallel for if (n_elts > CS_THR_MIN)
581 for (cs_lnum_t i = 0; i < n_elts; i++) {
582 const cs_lnum_t elt_id = elt_ids[i];
583 _icodcl[elt_id] = bc_type;
584 _rcodcl3[elt_id] = eval_buf[_dim*i + coo_id];
585 }
586
587 }
588 }
589 break;
590 }
591 }
592
593 /*============================================================================
594 * Fortran wrapper function definitions
595 *============================================================================*/
596
597 /*----------------------------------------------------------------------------
598 * Set mapped boundary conditions for a given field and mapping locator.
599 *
600 * parameters:
601 * field_id <-- id of field whose boundary conditions are set
602 * locator <-- associated mapping locator, as returned
603 * by cs_boundary_conditions_map().
604 * location_type <-- matching values location (CS_MESH_LOCATION_CELLS or
605 * CS_MESH_LOCATION_BOUNDARY_FACES)
606 * enforce_balance <-- balance handling option:
607 * 0: values are simply mapped
608 * 1: values are mapped, then multiplied
609 * by a constant factor so that their
610 * surface integral on selected faces
611 * is preserved (relative to the
612 * input values)
613 * 2: as 1, but with a boundary-defined
614 * weight, defined by balance_w
615 * 3: as 1, but with a cell-defined
616 * weight, defined by balance_w
617 * interpolate <-- interpolation option:
618 * 0: values are simply based on matching
619 * cell or face center values
620 * 1: values are based on matching cell
621 * or face center values, corrected
622 * by gradient interpolation
623 * n_faces <-- number of selected boundary faces
624 * faces <-- list of selected boundary faces (1 to n),
625 * or NULL if no indirection is needed
626 * balance_w <-- optional balance weight, or NULL
627 * nvar <-- number of variables requiring BC's
628 * rcodcl <-> boundary condition values
629 *----------------------------------------------------------------------------*/
630
631 void
cs_f_boundary_conditions_mapped_set(int field_id,ple_locator_t * locator,cs_mesh_location_type_t location_type,int enforce_balance,int interpolate,cs_lnum_t n_faces,const cs_lnum_t * faces,cs_real_t * balance_w,int nvar,cs_real_t * rcodcl)632 cs_f_boundary_conditions_mapped_set(int field_id,
633 ple_locator_t *locator,
634 cs_mesh_location_type_t location_type,
635 int enforce_balance,
636 int interpolate,
637 cs_lnum_t n_faces,
638 const cs_lnum_t *faces,
639 cs_real_t *balance_w,
640 int nvar,
641 cs_real_t *rcodcl)
642 {
643 cs_lnum_t *_faces = NULL;
644
645 if (faces != NULL) {
646 BFT_MALLOC(_faces, n_faces, cs_lnum_t);
647 for (cs_lnum_t i = 0; i < n_faces; i++)
648 _faces[i] = faces[i] - 1;
649 }
650
651 cs_boundary_conditions_mapped_set(cs_field_by_id(field_id),
652 locator,
653 location_type,
654 enforce_balance,
655 interpolate,
656 n_faces,
657 _faces,
658 balance_w,
659 nvar,
660 rcodcl);
661
662 BFT_FREE(_faces);
663 }
664
665 /*----------------------------------------------------------------------------
666 * Get pointer to cs_glob_bc_type and cs_glob_bc_face_zone
667 *----------------------------------------------------------------------------*/
668
669 void
cs_f_boundary_conditions_get_pointers(int ** itypfb,int ** izfppp)670 cs_f_boundary_conditions_get_pointers(int **itypfb,
671 int **izfppp)
672 {
673 *itypfb = _bc_type;
674 *izfppp = _bc_face_zone;
675 }
676
677 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
678
679 /*============================================================================
680 * Public function definitions
681 *============================================================================*/
682
683 /*----------------------------------------------------------------------------*/
684 /*!
685 * \brief Handling of boundary condition definition errors and
686 * associated output.
687 *
688 * This function checks for errors, and simply returns if no error is
689 * encountered. In case of error, it outputs helpful information so as to
690 * make it easier to locate the matching faces.
691 *
692 * For each boundary face, bc_type defines the boundary condition type.
693 * As a convention here, zero values correspond to undefined types,
694 * positive values to defined types (with no error), and negative values
695 * to defined types with inconsistent or incompatible values, the
696 * absolute value indicating the original boundary condition type.
697 *
698 * An optional label may be used if the error is related to another
699 * attribute than the boundary type, for appropriate error reporting.
700 *
701 * \param[in] bc_flag array of BC type ids
702 * \param[in] type_name name of attribute in error, or NULL
703 */
704 /*----------------------------------------------------------------------------*/
705
706 void
cs_boundary_conditions_error(const int * bc_flag,const char * type_name)707 cs_boundary_conditions_error(const int *bc_flag,
708 const char *type_name)
709 {
710 const char type_name_default[] = N_("boundary condition type");
711 const char *_type_name = type_name_default;
712
713 if (type_name != NULL)
714 _type_name = type_name;
715
716 /* Check for faces with problems;
717 bc_flag[] used to determine if there is an error */
718
719 int have_errors
720 = cs_flag_check(_("face with boundary condition definition error"),
721 _type_name,
722 _("BC type"),
723 _("Faces with B.C. error"),
724 _("Faces with valid B.C.'s"),
725 CS_MESH_LOCATION_BOUNDARY_FACES,
726 1, /* min_flag */
727 bc_flag);
728
729 if (have_errors)
730 bft_error
731 (__FILE__, __LINE__, 0,
732 _("\nSome boundary condition definitions are incomplete or incorrect.\n"
733 "\n For details, read the end of the calculation log,\n"
734 " or visualize the error postprocessing output."));
735 }
736
737 /*----------------------------------------------------------------------------*/
738 /*!
739 * \brief Locate shifted boundary face coordinates on possibly filtered
740 * cells or boundary faces for later interpolation.
741 *
742 * \param[in] location_type matching values location (CS_MESH_LOCATION_CELLS
743 * or CS_MESH_LOCATION_BOUNDARY_FACES)
744 * \param[in] n_location_elts number of selected location elements
745 * \param[in] n_faces number of selected boundary faces
746 * \param[in] location_elts list of selected location elements (0 to n-1),
747 * or NULL if no indirection is needed
748 * \param[in] faces list of selected boundary faces (0 to n-1),
749 * or NULL if no indirection is needed
750 * \param[in] coord_shift array of coordinates shift relative to selected
751 * boundary faces
752 * \param[in] coord_stride access stride in coord_shift: 0 for uniform
753 * shift, 1 for "per face" shift.
754 * \param[in] tolerance relative tolerance for point location.
755 *
756 * \return associated locator structure
757 */
758 /*----------------------------------------------------------------------------*/
759
760 ple_locator_t *
cs_boundary_conditions_map(cs_mesh_location_type_t location_type,cs_lnum_t n_location_elts,cs_lnum_t n_faces,const cs_lnum_t * location_elts,const cs_lnum_t * faces,cs_real_3_t * coord_shift,int coord_stride,double tolerance)761 cs_boundary_conditions_map(cs_mesh_location_type_t location_type,
762 cs_lnum_t n_location_elts,
763 cs_lnum_t n_faces,
764 const cs_lnum_t *location_elts,
765 const cs_lnum_t *faces,
766 cs_real_3_t *coord_shift,
767 int coord_stride,
768 double tolerance)
769 {
770 ple_locator_t *locator = NULL;
771
772 /* Build temporary "donor" location mesh */
773 /*----------------------------------------*/
774
775 fvm_nodal_t *nm = NULL;
776
777 /* Convert from 0-based to 1-based numbering (note that in practice,
778 we probably use 1-based numbering upstream, and convert back to
779 1-based downstream, but we prefer to use a 0-based API here so as
780 to make progress on global 1-based to 0-based conversion). */
781
782 cs_lnum_t *_location_elts = NULL;
783 if (location_elts != NULL) {
784 BFT_MALLOC(_location_elts, n_location_elts, cs_lnum_t);
785 for (cs_lnum_t i = 0; i < n_location_elts; i++)
786 _location_elts[i] = location_elts[i] + 1;
787 }
788
789 if (location_type == CS_MESH_LOCATION_CELLS)
790 nm = cs_mesh_connect_cells_to_nodal(cs_glob_mesh,
791 "search mesh",
792 false, /* include_families */
793 n_location_elts,
794 _location_elts);
795 else if (location_type == CS_MESH_LOCATION_BOUNDARY_FACES)
796 nm = cs_mesh_connect_faces_to_nodal(cs_glob_mesh,
797 "search mesh",
798 false, /* include_families */
799 0,
800 n_location_elts,
801 NULL,
802 _location_elts);
803
804 BFT_FREE(_location_elts);
805
806 /* Now build locator */
807 /*-------------------*/
808
809 #if defined(PLE_HAVE_MPI)
810 locator = ple_locator_create(cs_glob_mpi_comm, cs_glob_n_ranks, 0);
811 #else
812 locator = ple_locator_create();
813 #endif
814
815 int options[PLE_LOCATOR_N_OPTIONS];
816 for (int i = 0; i < PLE_LOCATOR_N_OPTIONS; i++)
817 options[i] = 0;
818 options[PLE_LOCATOR_NUMBERING] = 0; /* base 0 numbering */
819
820 /* Build location coordinates */
821
822 ple_coord_t *point_coords;
823
824 const cs_real_3_t *restrict b_face_cog
825 = (const cs_real_3_t *restrict)cs_glob_mesh_quantities->b_face_cog;
826
827 BFT_MALLOC(point_coords, n_faces*3, ple_coord_t);
828 if (faces != NULL) {
829 for (cs_lnum_t i = 0; i < n_faces; i++) {
830 const cs_lnum_t face_id = faces[i];
831 for (cs_lnum_t j = 0; j < 3; j++)
832 point_coords[i*3 + j] = b_face_cog[face_id][j]
833 + coord_shift[i*coord_stride][j];
834 }
835 }
836
837 ple_locator_set_mesh(locator,
838 nm,
839 options,
840 0., /* tolerance_base */
841 tolerance,
842 3, /* dim */
843 n_faces,
844 NULL,
845 NULL, /* point_tag */
846 point_coords,
847 NULL, /* distance */
848 cs_coupling_mesh_extents,
849 cs_coupling_point_in_mesh_p);
850
851 /* Check that location is OK */
852
853 cs_gnum_t loc_count[2];
854 loc_count[0] = ple_locator_get_n_exterior(locator);
855 loc_count[1] = ple_locator_get_n_exterior(locator);
856 cs_parall_counter(loc_count, 2);
857
858 if (loc_count[1] > 0) {
859 bft_error
860 (__FILE__, __LINE__, 0,
861 _("\nIn function %s,\n"
862 " %llu boundary faces (of %llu selected) were not matched to mesh\n"
863 " elements. Check your coordinate shift definitions."),
864 __func__,
865 (unsigned long long)loc_count[1],
866 (unsigned long long)loc_count[0]);
867 }
868
869 BFT_FREE(point_coords);
870
871 /* Shift from 1-base to 0-based locations */
872
873 ple_locator_shift_locations(locator, -1);
874
875 /* Nodal mesh is not needed anymore */
876
877 nm = fvm_nodal_destroy(nm);
878
879 /* Return initialized locator */
880
881 return locator;
882 }
883
884 /*----------------------------------------------------------------------------*/
885 /*!
886 * \brief Set mapped boundary conditions for a given field and mapping locator.
887 *
888 * \param[in] f field whose boundary conditions are set
889 * \param[in] locator associated mapping locator, as returned
890 * by \ref cs_boundary_conditions_map.
891 * \param[in] location_type matching values location
892 * (CS_MESH_LOCATION_CELLS or
893 * CS_MESH_LOCATION_BOUNDARY_FACES)
894 * \param[in] normalize normalization option:
895 * 0: values are simply mapped
896 * 1: values are mapped, then multiplied
897 * by a constant factor so that their
898 * surface integral on selected faces
899 * is preserved (relative to the
900 * input values)
901 * 2: as 1, but with a boundary-defined
902 * weight, defined by balance_w
903 * 3: as 1, but with a cell-defined
904 * weight, defined by balance_w
905 * \param[in] interpolate interpolation option:
906 * 0: values are simply based on matching
907 * cell or face center values
908 * 1: values are based on matching cell
909 * or face center values, corrected
910 * by gradient interpolation
911 * \param[in] n_faces number of selected boundary faces
912 * \param[in] faces list of selected boundary faces (0 to n-1),
913 * or NULL if no indirection is needed
914 * \param[in] balance_w optional balance weight, or NULL
915 * \param[in] nvar number of variables requiring BC's
916 * \param[in, out] rcodcl boundary condition values
917 */
918 /*----------------------------------------------------------------------------*/
919
920 void
cs_boundary_conditions_mapped_set(const cs_field_t * f,ple_locator_t * locator,cs_mesh_location_type_t location_type,int normalize,int interpolate,cs_lnum_t n_faces,const cs_lnum_t * faces,cs_real_t * balance_w,int nvar,cs_real_t rcodcl[])921 cs_boundary_conditions_mapped_set(const cs_field_t *f,
922 ple_locator_t *locator,
923 cs_mesh_location_type_t location_type,
924 int normalize,
925 int interpolate,
926 cs_lnum_t n_faces,
927 const cs_lnum_t *faces,
928 cs_real_t *balance_w,
929 int nvar,
930 cs_real_t rcodcl[])
931 {
932 int var_id = -1;
933
934 const int dim = f->dim;
935 const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
936 const ple_lnum_t n_dist = ple_locator_get_n_dist_points(locator);
937 const ple_lnum_t *dist_loc = ple_locator_get_dist_locations(locator);
938 const ple_coord_t *dist_coords = ple_locator_get_dist_coords(locator);
939
940 cs_field_interpolate_t interpolation_type = CS_FIELD_INTERPOLATE_MEAN;
941
942 cs_real_t inlet_sum_0[9], inlet_sum_1[9];
943 cs_real_t *distant_var, *local_var;
944
945 /* Get field's variable id */
946
947 static int var_id_key = -1;
948 if (var_id_key < 0)
949 var_id_key = cs_field_key_id("variable_id");
950 assert(var_id_key >= 0);
951
952 var_id = cs_field_get_key_int(f, var_id_key) - 1;
953
954 if (var_id < 0)
955 return;
956
957 assert(f->location_id == CS_MESH_LOCATION_CELLS);
958
959 /* Compute initial balance if applicable */
960
961 if (normalize > 0) {
962 assert(dim <= 9);
963 _inlet_sum(var_id,
964 f,
965 cs_glob_mesh,
966 cs_glob_mesh_quantities,
967 normalize,
968 n_faces,
969 faces,
970 balance_w,
971 nvar,
972 rcodcl,
973 inlet_sum_0);
974 }
975
976 /* Allocate working array */
977
978 BFT_MALLOC(distant_var, n_dist*dim, cs_real_t);
979 BFT_MALLOC(local_var, n_faces*dim, cs_real_t);
980
981 /* Prepare values to send */
982 /*------------------------*/
983
984 if (interpolate)
985 interpolation_type = CS_FIELD_INTERPOLATE_GRADIENT;
986
987 assert(sizeof(ple_coord_t) == sizeof(cs_real_t));
988
989 if (location_type == CS_MESH_LOCATION_CELLS || interpolate) {
990 /* FIXME: we cheat here with the constedness of the field
991 for a possible ghost values update, but having a finer control
992 of when syncing is required would be preferable */
993 cs_field_t *_f = cs_field_by_id(f->id);
994 cs_field_interpolate(_f,
995 interpolation_type,
996 n_dist,
997 dist_loc,
998 (const cs_real_3_t *)dist_coords,
999 distant_var);
1000 }
1001
1002 else if (location_type == CS_MESH_LOCATION_BOUNDARY_FACES) {
1003
1004 const cs_lnum_t *b_face_cells = cs_glob_mesh->b_face_cells;
1005 const cs_field_bc_coeffs_t *bc_coeffs = f->bc_coeffs;
1006
1007 /* If boundary condition coefficients are available */
1008
1009 if (bc_coeffs != NULL) {
1010
1011 if (dim == 1) {
1012 for (cs_lnum_t i = 0; i < n_dist; i++) {
1013 cs_lnum_t f_id = dist_loc[i];
1014 cs_lnum_t c_id = b_face_cells[f_id];
1015 distant_var[i] = bc_coeffs->a[f_id]
1016 + bc_coeffs->b[f_id]*f->val[c_id];
1017 }
1018 }
1019 else {
1020 for (cs_lnum_t i = 0; i < n_dist; i++) {
1021 cs_lnum_t f_id = dist_loc[i];
1022 cs_lnum_t c_id = b_face_cells[f_id];
1023 for (cs_lnum_t j = 0; j < dim; j++) {
1024 distant_var[i*dim+j] = bc_coeffs->a[f_id*dim+j];
1025 for (cs_lnum_t k = 0; k < dim; k++)
1026 distant_var[i*dim+j] += bc_coeffs->b[(f_id*dim+k)*dim + j]
1027 *f->val[c_id*dim+k];
1028 }
1029 }
1030 }
1031
1032 }
1033
1034 /* If no boundary condition coefficients are available */
1035
1036 else {
1037
1038 for (cs_lnum_t i = 0; i < n_dist; i++) {
1039 cs_lnum_t f_id = dist_loc[i];
1040 cs_lnum_t c_id = b_face_cells[f_id];
1041 for (cs_lnum_t j = 0; j < dim; j++)
1042 distant_var[i*dim+j] = f->val[c_id*dim+j];
1043 }
1044
1045 }
1046
1047 }
1048
1049 ple_locator_exchange_point_var(locator,
1050 distant_var,
1051 local_var,
1052 NULL, /* faces indirection */
1053 sizeof(cs_real_t),
1054 f->dim,
1055 0);
1056
1057 /* Now set boundary condition values */
1058
1059 for (cs_lnum_t j = 0; j < dim; j++) {
1060
1061 cs_real_t *_rcodcl = rcodcl + (var_id+j)*n_b_faces;
1062
1063 for (cs_lnum_t i = 0; i < n_faces; i++) {
1064 const cs_lnum_t f_id = (faces != NULL) ? faces[i] : i;
1065 _rcodcl[f_id] = local_var[i*dim + j];
1066 }
1067
1068 }
1069
1070 BFT_FREE(local_var);
1071 BFT_FREE(distant_var);
1072
1073 /* Compute initial balance if applicable */
1074
1075 if (normalize > 0) {
1076
1077 _inlet_sum(var_id,
1078 f,
1079 cs_glob_mesh,
1080 cs_glob_mesh_quantities,
1081 normalize,
1082 n_faces,
1083 faces,
1084 balance_w,
1085 nvar,
1086 rcodcl,
1087 inlet_sum_1);
1088
1089 for (cs_lnum_t j = 0; j < dim; j++) {
1090
1091 const cs_real_t f_mult = (fabs(inlet_sum_1[j]) > 1.e-24) ?
1092 inlet_sum_0[j] / inlet_sum_1[j] : 1.;
1093
1094 cs_real_t *_rcodcl = rcodcl + (var_id+j)*n_b_faces;
1095
1096 for (cs_lnum_t i = 0; i < n_faces; i++) {
1097 const cs_lnum_t f_id = (faces != NULL) ? faces[i] : i;
1098 _rcodcl[f_id] *= f_mult;
1099 }
1100
1101 }
1102
1103 }
1104 }
1105
1106 /*----------------------------------------------------------------------------*/
1107 /*!
1108 * \brief Add location of locate shifted boundary face coordinates on
1109 * cells or boundary faces for automatic interpolation.
1110 *
1111 * \note
1112 * This function is currently restricted to mapping of boundary face
1113 * locations (usually from boundary zones) to cell of boundary face
1114 * locations, but could be extended to other location types in the future.
1115 *
1116 * \param[in] bc_location_id id of selected boundary mesh location;
1117 * currently restricted to subsets of
1118 * boundary faces (i.e. boundary zone
1119 * location ids).
1120 * \param[in] source_location_id id of selected location mesh location
1121 * (usually CS_MESH_LOCATION_CELLS but can be
1122 * a more restricted cell or boundary face zone
1123 * location location id).
1124 * \param[in] coord_shift coordinates shift relative to selected
1125 * boundary faces
1126 * \param[in] tolerance relative tolerance for point location.
1127 *
1128 * \return id of added map
1129 */
1130 /*----------------------------------------------------------------------------*/
1131
1132 int
cs_boundary_conditions_add_map(int bc_location_id,int source_location_id,cs_real_t coord_shift[3],double tolerance)1133 cs_boundary_conditions_add_map(int bc_location_id,
1134 int source_location_id,
1135 cs_real_t coord_shift[3],
1136 double tolerance)
1137 {
1138 int map_id = _n_bc_maps;
1139
1140 BFT_REALLOC(_bc_maps, _n_bc_maps+1, cs_bc_map_t);
1141
1142 cs_bc_map_t *bc_map = _bc_maps + _n_bc_maps;
1143
1144 bc_map->bc_location_id = bc_location_id;
1145 bc_map->source_location_id = source_location_id;
1146
1147 for (int i = 0; i < 3; i++)
1148 bc_map->coord_shift[i] = coord_shift[i];
1149
1150 bc_map->tolerance = tolerance;
1151
1152 bc_map->locator = NULL;
1153
1154 _n_bc_maps += 1;
1155
1156 return map_id;
1157 }
1158
1159 /*----------------------------------------------------------------------------*/
1160 /*!
1161 * \brief Create the legacy boundary conditions face type and face zone arrays.
1162 */
1163 /*----------------------------------------------------------------------------*/
1164
1165 void
cs_boundary_conditions_create(void)1166 cs_boundary_conditions_create(void)
1167 {
1168 const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
1169
1170 /* Get default boundary type (converting "current" to "legacy" codes) */
1171
1172 const cs_boundary_t *boundaries = cs_glob_boundaries;
1173 int default_type = 0;
1174 if (boundaries->default_type & CS_BOUNDARY_WALL)
1175 default_type = CS_SMOOTHWALL;
1176 else if (boundaries->default_type & CS_BOUNDARY_SYMMETRY)
1177 default_type = CS_SYMMETRY;
1178
1179 /* boundary conditions type by boundary face */
1180
1181 BFT_MALLOC(_bc_type, n_b_faces, int);
1182 for (cs_lnum_t ii = 0 ; ii < n_b_faces ; ii++) {
1183 _bc_type[ii] = default_type;
1184 }
1185 cs_glob_bc_type = _bc_type;
1186
1187 /* boundary conditions zone by boundary face */
1188 /* only for specific physics */
1189
1190 BFT_MALLOC(_bc_face_zone, n_b_faces, int);
1191 for (cs_lnum_t ii = 0 ; ii < n_b_faces ; ii++) {
1192 _bc_face_zone[ii] = 0;
1193 }
1194 cs_glob_bc_face_zone = _bc_face_zone;
1195 }
1196
1197 /*----------------------------------------------------------------------------*/
1198 /*!
1199 * \brief Free the boundary conditions face type and face zone arrays.
1200 *
1201 * This also frees boundary condition mappings which may have been defined.
1202 */
1203 /*----------------------------------------------------------------------------*/
1204
1205 void
cs_boundary_conditions_free(void)1206 cs_boundary_conditions_free(void)
1207 {
1208 BFT_FREE(_bc_type);
1209 BFT_FREE(_bc_face_zone);
1210
1211 for (int i = 0; i < _n_bc_maps; i++)
1212 ple_locator_destroy((_bc_maps + i)->locator);
1213
1214 BFT_FREE(_bc_maps);
1215 _n_bc_maps = 0;
1216 }
1217
1218 /*----------------------------------------------------------------------------*/
1219 /*!
1220 * \brief Set convective oulet boundary condition for a scalar.
1221 *
1222 * \param[out] coefa explicit BC coefficient for gradients
1223 * \param[out] cofaf explicit BC coefficient for diffusive flux
1224 * \param[out] coefb implicit BC coefficient for gradients
1225 * \param[out] cofbf implicit BC coefficient for diffusive flux
1226 * \param[in] pimp Flux value to impose
1227 * \param[in] cfl Local Courant number used to convect
1228 * \param[in] hint Internal exchange coefficient
1229 */
1230 /*----------------------------------------------------------------------------*/
1231
1232 void
cs_boundary_conditions_set_convective_outlet_scalar(cs_real_t * coefa,cs_real_t * cofaf,cs_real_t * coefb,cs_real_t * cofbf,cs_real_t pimp,cs_real_t cfl,cs_real_t hint)1233 cs_boundary_conditions_set_convective_outlet_scalar(cs_real_t *coefa ,
1234 cs_real_t *cofaf,
1235 cs_real_t *coefb,
1236 cs_real_t *cofbf,
1237 cs_real_t pimp,
1238 cs_real_t cfl,
1239 cs_real_t hint)
1240 {
1241 /* Gradient BCs */
1242 *coefb = cfl / (1.0 + cfl);
1243 *coefa = (1.0 - *coefb) * pimp;
1244
1245 /* Flux BCs */
1246 *cofaf = - hint * *coefa;
1247 *cofbf = hint * (1.0 - *coefb);
1248 }
1249
1250 /*----------------------------------------------------------------------------*/
1251 /*!
1252 * \brief Update per variable boundary condition codes.
1253 *
1254 * \param[in] nvar number of variables requiring BC's
1255 * \param[in] itypfb type of boundary for each face
1256 * \param[in, out] icodcl boundary condition codes
1257 * \param[in, out] rcodcl boundary condition values
1258 */
1259 /*----------------------------------------------------------------------------*/
1260
1261 void
cs_boundary_conditions_compute(int nvar,int itypfb[],int icodcl[],double rcodcl[])1262 cs_boundary_conditions_compute(int nvar,
1263 int itypfb[],
1264 int icodcl[],
1265 double rcodcl[])
1266 {
1267 CS_UNUSED(itypfb);
1268
1269 /* Initialization */
1270
1271 const cs_time_step_t *ts = cs_glob_time_step;
1272 const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
1273
1274 for (int map_id = 0; map_id < _n_bc_maps; map_id++)
1275 _update_bc_map(map_id);
1276
1277 static int var_id_key = -1;
1278 if (var_id_key < 0)
1279 var_id_key = cs_field_key_id("variable_id");
1280 assert(var_id_key >= 0);
1281
1282 const cs_mesh_t *mesh = cs_glob_mesh;
1283 const cs_boundary_t *boundaries = cs_glob_boundaries;
1284 const cs_real_t t_eval = ts->t_cur;
1285
1286 /* Buffer for evaluation of analytical values */
1287
1288 cs_lnum_t eval_buf_size = 0;
1289 cs_real_t *eval_buf = NULL;
1290
1291 /* Loop on fields */
1292
1293 const int n_fields = cs_field_n_fields();
1294
1295 for (int f_id = 0; f_id < n_fields; f_id++) {
1296 cs_field_t *f = cs_field_by_id(f_id);
1297 const cs_equation_param_t *eqp = NULL;
1298 if (f->type & CS_FIELD_VARIABLE)
1299 eqp = cs_field_get_equation_param_const(f);
1300
1301 if (eqp == NULL)
1302 continue;
1303
1304 /* Only handle legacy discretization here */
1305
1306 if (eqp->space_scheme != CS_SPACE_SCHEME_LEGACY)
1307 continue;
1308
1309 /* Settings in eqp may not be well-defined at this stage. The following
1310 test should be more appropriate to decide if one skips this field or
1311 not */
1312 if (f->type & CS_FIELD_CDO)
1313 continue;
1314
1315 char description[128];
1316 snprintf(description, 127, _("boundary condition definitions for \"%s\""),
1317 f->name);
1318 description[127] = '\0';
1319
1320 /* Get associated variable id */
1321
1322 int var_id = cs_field_get_key_int(f, var_id_key) - 1;
1323 assert(var_id >= 0);
1324
1325 /* Conversion flag for enthalpy (temperature given);
1326 not active yet in GUI and XML. */
1327
1328 int icodcl_m = 1;
1329
1330 cs_lnum_t n_max_vals = (cs_lnum_t)(f->dim) * n_b_faces;
1331
1332 if (n_max_vals > eval_buf_size) {
1333 eval_buf_size = n_max_vals;
1334 BFT_FREE(eval_buf);
1335 BFT_MALLOC(eval_buf, eval_buf_size, cs_real_t);
1336 }
1337
1338 /* Loop on boundary conditions */
1339
1340 for (int bc_id = 0; bc_id < eqp->n_bc_defs; bc_id++) {
1341 const cs_xdef_t *def = eqp->bc_defs[bc_id];
1342 cs_param_bc_type_t bc_type = (cs_param_bc_type_t)(def->meta);
1343 switch (bc_type) {
1344
1345 case CS_PARAM_BC_HMG_DIRICHLET:
1346 _compute_hmg_dirichlet_bc(mesh,
1347 boundaries,
1348 eqp,
1349 def,
1350 var_id,
1351 icodcl_m,
1352 icodcl,
1353 rcodcl);
1354 break;
1355
1356 case CS_PARAM_BC_DIRICHLET:
1357 _compute_dirichlet_bc(mesh,
1358 boundaries,
1359 f,
1360 eqp,
1361 def,
1362 description,
1363 var_id,
1364 t_eval,
1365 icodcl_m,
1366 icodcl,
1367 rcodcl,
1368 eval_buf);
1369 break;
1370
1371 case CS_PARAM_BC_HMG_NEUMANN:
1372 _compute_hmg_neumann_bc(mesh,
1373 def,
1374 nvar,
1375 var_id,
1376 icodcl,
1377 rcodcl);
1378 break;
1379
1380 case CS_PARAM_BC_NEUMANN:
1381 _compute_neumann_bc(mesh,
1382 eqp,
1383 def,
1384 description,
1385 nvar,
1386 var_id,
1387 t_eval,
1388 icodcl,
1389 rcodcl,
1390 eval_buf);
1391 break;
1392
1393 default:
1394 break;
1395 }
1396 }
1397
1398 }
1399
1400 BFT_FREE(eval_buf);
1401 }
1402
1403 /*----------------------------------------------------------------------------*/
1404 /*!
1405 * \brief Automatic adjustments for boundary condition codes.
1406 *
1407 * Currently handles mapped inlets, after the call to \ref stdtcl.
1408 * As portions of stdtcl are migrated to C, they should be called here,
1409 * before mapped inlets.
1410 *
1411 * \param[in] nvar number of variables requiring BC's
1412 * \param[in] itypfb type of boundary for each face
1413 * \param[in, out] icodcl boundary condition codes
1414 * \param[in, out] rcodcl boundary condition values
1415 */
1416 /*----------------------------------------------------------------------------*/
1417
1418 void
cs_boundary_conditions_complete(int nvar,int itypfb[],int icodcl[],double rcodcl[])1419 cs_boundary_conditions_complete(int nvar,
1420 int itypfb[],
1421 int icodcl[],
1422 double rcodcl[])
1423 {
1424 CS_NO_WARN_IF_UNUSED(itypfb);
1425
1426 /* Initialization */
1427
1428 const cs_time_step_t *ts = cs_glob_time_step;
1429
1430 for (int map_id = 0; map_id < _n_bc_maps; map_id++)
1431 _update_bc_map(map_id);
1432
1433 static int var_id_key = -1;
1434 if (var_id_key < 0)
1435 var_id_key = cs_field_key_id("variable_id");
1436 assert(var_id_key >= 0);
1437
1438 /* Loop on fields */
1439
1440 const int n_fields = cs_field_n_fields();
1441
1442 for (int f_id = 0; f_id < n_fields; f_id++) {
1443 cs_field_t *f = cs_field_by_id(f_id);
1444
1445 if (! (f->type & CS_FIELD_VARIABLE))
1446 continue;
1447
1448 const cs_equation_param_t *eqp
1449 = cs_field_get_equation_param_const(f);
1450
1451 /* Only handle legacy discretization here */
1452 if (eqp != NULL) {
1453 if (eqp->space_scheme != CS_SPACE_SCHEME_LEGACY)
1454 continue;
1455 }
1456
1457 /* Settings in eqp may not be well-defined at this stage. The following
1458 test should be more appropriate to decide if one skips this field or
1459 not */
1460 if (f->type & CS_FIELD_CDO)
1461 continue;
1462
1463 /* Get associated variable id */
1464
1465 int var_id = cs_field_get_key_int(f, var_id_key) - 1;
1466 assert(var_id >= 0);
1467
1468 /* Treatment of mapped inlets */
1469
1470 for (int map_id = 0; map_id < _n_bc_maps; map_id++) {
1471
1472 cs_bc_map_t *bc_map = _bc_maps + map_id;
1473
1474 if (bc_map->locator == NULL || ts->nt_cur <= 1)
1475 continue;
1476
1477 int interpolate = 0;
1478 int normalize = 0;
1479 if (f == CS_F_(vel))
1480 normalize = 1;
1481 else {
1482 const int keysca = cs_field_key_id("scalar_id");
1483 if (cs_field_get_key_int(f, keysca) > 0)
1484 normalize = 1;
1485 }
1486
1487 if (f != CS_F_(p)) {
1488 cs_mesh_location_type_t location_type
1489 = cs_mesh_location_get_type(bc_map->source_location_id);
1490 cs_lnum_t n_faces
1491 = cs_mesh_location_get_n_elts(bc_map->bc_location_id)[0];
1492 const cs_lnum_t *faces
1493 = cs_mesh_location_get_elt_ids_try(bc_map->bc_location_id);
1494
1495 cs_boundary_conditions_mapped_set(f,
1496 bc_map->locator,
1497 location_type,
1498 normalize,
1499 interpolate,
1500 n_faces,
1501 faces,
1502 NULL,
1503 nvar,
1504 rcodcl);
1505
1506 if (f == CS_F_(h)) {
1507 const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
1508 int *_icodcl = icodcl + var_id*n_b_faces;
1509
1510 for (cs_lnum_t i = 0; i < n_faces; i++) {
1511 const cs_lnum_t j = (faces != NULL) ? faces[i] : i;
1512 if (_icodcl[j] < 0)
1513 _icodcl[j] *= -1;
1514 }
1515 }
1516
1517 }
1518
1519 }
1520
1521 }
1522 }
1523
1524 /*----------------------------------------------------------------------------*/
1525
1526 END_C_DECLS
1527