1 /*============================================================================
2 * Management of the GUI parameters file: boundary conditions
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <stdio.h>
34 #include <stdlib.h>
35 #include <math.h>
36 #include <string.h>
37 #include <fcntl.h>
38 #include <unistd.h>
39 #include <assert.h>
40
41 /*----------------------------------------------------------------------------
42 * Local headers
43 *----------------------------------------------------------------------------*/
44
45 #include "bft_mem.h"
46 #include "bft_error.h"
47 #include "bft_printf.h"
48
49 #include "fvm_selector.h"
50
51 #include "cs_ale.h"
52 #include "cs_atmo.h"
53 #include "cs_base.h"
54 #include "cs_boundary.h"
55 #include "cs_boundary_conditions.h"
56 #include "cs_boundary_zone.h"
57 #include "cs_combustion_model.h"
58 #include "cs_equation_param.h"
59 #include "cs_parameters.h"
60 #include "cs_gui_util.h"
61 #include "cs_gui.h"
62 #include "cs_gui_specific_physics.h"
63 #include "cs_mesh.h"
64 #include "cs_field.h"
65 #include "cs_field_default.h"
66 #include "cs_field_pointer.h"
67 #include "cs_physical_model.h"
68 #include "cs_thermal_model.h"
69 #include "cs_timer.h"
70 #include "cs_tree.h"
71 #include "cs_turbulence_model.h"
72 #include "cs_parall.h"
73 #include "cs_elec_model.h"
74 #include "cs_physical_model.h"
75 #include "cs_prototypes.h"
76 #include "cs_wall_functions.h"
77
78 /*----------------------------------------------------------------------------
79 * Header for the current file
80 *----------------------------------------------------------------------------*/
81
82 #include "cs_gui_boundary_conditions.h"
83
84 /*----------------------------------------------------------------------------*/
85
86 BEGIN_C_DECLS
87
88 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
89
90 /*=============================================================================
91 * Local Macro Definitions
92 *============================================================================*/
93
94 /* debugging switch */
95 #define _XML_DEBUG_ 0
96
97 /*=============================================================================
98 * Local Type Definitions
99 *============================================================================*/
100
101 /* Enum for boundary conditions */
102
103 typedef enum {
104 BY_XDEF = -1, /* to mark usage of newer system */
105 DIRICHLET_CNV,
106 DIRICHLET_FORMULA,
107 DIRICHLET_IMPLICIT,
108 EXCHANGE_COEFF,
109 EXCHANGE_COEFF_FORMULA,
110 FLOW1,
111 HYDRAULIC_DIAMETER,
112 NEUMANN_FORMULA,
113 NEUMANN_IMPLICIT,
114 TURBULENT_INTENSITY,
115 } cs_boundary_value_t;
116
117 typedef struct {
118 double val1; /* fortran array RCODCL(.,.,1) mapping */
119 double val2; /* fortran array RCODCL(.,.,2) mapping */
120 } cs_val_t;
121
122 typedef struct {
123 int read_data; /* 1 if profile is calculated from data */
124 int automatic; /* 1 if nature of the boundary is automatic */
125 } cs_meteo_t;
126
127 typedef struct {
128
129 int n_fields; /* number of handled fields */
130 int n_zones; /* number of associated zones */
131
132 const char **label; /* pointer to label for each boundary zone */
133 const char **nature; /* nature for each boundary zone */
134 int *bc_num; /* associated number */
135
136 int *iqimp; /* 1 if a flow rate is applied */
137 int *ientfu; /* 1 for a fuel flow inlet (gas combustion - D3P) */
138 int *ientox; /* 1 for an air flow inlet (gas combustion - D3P) */
139 int *ientgb; /* 1 for burned gas inlet (gas combustion) */
140 int *ientgf; /* 1 for unburned gas inlet (gas combustion) */
141 int *ientat; /* 1 if inlet for oxydant (coal combustion) */
142 int *ientcp; /* 1 if inlet for oxydant+coal (coal combustion) */
143 int *icalke; /* automatic boundaries for turbulent variables */
144 double *qimp; /* oxydant flow rate (coal combustion) */
145 int *inmoxy; /* oxydant number (coal combustion) */
146 double *timpat; /* inlet temperature of oxydant (coal combustion) */
147 double *tkent; /* inlet temperature (gas combustion) */
148 double **qimpcp; /* inlet coal flow rate (coal combustion) */
149 double **timpcp; /* inlet coal temperature (coal combustion) */
150 double *fment; /* Mean Mixture Fraction at Inlet (gas combustion) */
151 int *itype; /* type of inlet/outlet (compressible model) */
152 double *prein; /* inlet pressure (compressible model) */
153 double *rhoin; /* inlet density (compressible model) */
154 double *tempin; /* inlet temperature (compressible model) */
155 double *dh; /* inlet hydraulic diameter */
156 double *xintur; /* inlet turbulent intensity */
157 int **type_code; /* type of boundary for each variable */
158 cs_val_t **values; /* fortran array RCODCL mapping */
159 double ***distch; /* ratio for each coal */
160 double *rough; /* roughness size */
161 double *norm; /* norm of velocity vector */
162 cs_real_3_t *dir; /* directions inlet velocity */
163 bool *t_to_h; /* convert Enthalpy to temperature */
164 bool *velocity_e; /* formula for norm or mass flow rate of velocity */
165 bool *direction_e; /* formula for direction of velocity */
166 bool **scalar_e; /* formula for scalar (neumann, dirichlet or
167 exchange coefficient) */
168 bool *head_loss_e; /* formula for head loss (free inlet/outlet) */
169 bool *groundwat_e; /* formula for hydraulic head (groundwater) */
170
171 cs_meteo_t *meteo; /* inlet or outlet info for atmospheric flow */
172
173 } cs_gui_boundary_t;
174
175 /*============================================================================
176 * Static global variables
177 *============================================================================*/
178
179 /* Pointer on the main boundaries structure */
180
181 static cs_gui_boundary_t *boundaries = NULL;
182
183 /*============================================================================
184 * Private function definitions
185 *============================================================================*/
186
187 /*----------------------------------------------------------------------------
188 * Return a pointer to equation parameters based on a field or equation name.
189 *
190 * parameters:
191 * name <-- field or equation name
192 *
193 * return:
194 * pointer to matching child string
195 *----------------------------------------------------------------------------*/
196
197 static cs_equation_param_t *
_get_equation_param(const char * name)198 _get_equation_param(const char *name)
199 {
200 cs_equation_param_t *eqp = NULL;
201
202 cs_field_t *f = cs_field_by_name_try(name);
203 if (f != NULL)
204 eqp = cs_field_get_equation_param(f);
205
206 /* FIXME: else get by equation name */
207
208 return eqp;
209 }
210
211 /*----------------------------------------------------------------------------
212 * Return a node associated with a given zone's boundary condition definition.
213 *
214 * parameters:
215 * tn <-- first node associated with search
216 *
217 * return:
218 * pointer to matching child string
219 *----------------------------------------------------------------------------*/
220
221 static cs_tree_node_t *
_get_zone_bc_node(cs_tree_node_t * tn_start,int izone)222 _get_zone_bc_node(cs_tree_node_t *tn_start,
223 int izone)
224 {
225 cs_tree_node_t *tn = tn_start;
226
227 /* if the start BC node is of a different type, search from parent */
228
229 if (strcmp(tn->name, boundaries->nature[izone]))
230 tn = cs_tree_node_get_child(tn_start->parent,
231 boundaries->nature[izone]);
232
233 /* Now searh from siblings */
234
235 tn = cs_tree_node_get_sibling_with_tag
236 (tn, "label", boundaries->label[izone]);
237
238 return tn;
239 }
240
241 /*-----------------------------------------------------------------------------
242 * Get status of data for inlet or outlet information.
243 *
244 * parameters:
245 * nature <-- nature of the boundary
246 * label <-- label of the inlet or the outlet
247 * tag <-- name of researched data
248 * data --> value associated to the data
249 *----------------------------------------------------------------------------*/
250
251 static void
_boundary_status(const char * nature,const char * label,const char * tag,int * data)252 _boundary_status(const char *nature,
253 const char *label,
254 const char *tag,
255 int *data)
256 {
257 cs_tree_node_t *tn
258 = cs_tree_get_node(cs_glob_tree, "boundary_conditions");
259 tn = cs_tree_get_node(tn, nature);
260 tn = cs_tree_node_get_sibling_with_tag(tn, "label", label);
261 tn = cs_tree_get_node(tn, tag);
262
263 cs_gui_node_get_status_int(tn, data);
264 }
265
266 /*-----------------------------------------------------------------------------
267 * Get status of data for velocity_pressure node of inlet or outlet information.
268 *
269 * parameters:
270 * nature <-- nature of the boundary
271 * label <-- label of the inlet or the outlet
272 * tag <-- name of researched data
273 * data --> value associated to the data
274 *----------------------------------------------------------------------------*/
275
276 static void
_boundary_status_vp(const char * nature,const char * label,const char * tag,int * data)277 _boundary_status_vp(const char *nature,
278 const char *label,
279 const char *tag,
280 int *data)
281 {
282 cs_tree_node_t *tn
283 = cs_tree_get_node(cs_glob_tree, "boundary_conditions");
284 tn = cs_tree_get_node(tn, nature);
285 tn = cs_tree_node_get_sibling_with_tag(tn, "label", label);
286 tn = cs_tree_get_node(tn, "velocity_pressure");
287 tn = cs_tree_get_node(tn, tag);
288 cs_gui_node_get_status_int(tn, data);
289 }
290
291 /*-----------------------------------------------------------------------------
292 * Check if a zone uses a mapped inlet, and define the associated mapping
293 * if this is the case.
294 *
295 * parameters:
296 * label <-- label of wall boundary condition
297 * z <-- pointer to boundary zone
298 *----------------------------------------------------------------------------*/
299
300 static void
_check_and_add_mapped_inlet(const char * label,const cs_zone_t * z)301 _check_and_add_mapped_inlet(const char *label,
302 const cs_zone_t *z)
303 {
304 int mapped_inlet = 0;
305
306 cs_tree_node_t *tn
307 = cs_tree_get_node(cs_glob_tree, "boundary_conditions");
308 tn = cs_tree_get_node(tn, "inlet");
309 tn = cs_tree_node_get_sibling_with_tag(tn, "label", label);
310
311 tn = cs_tree_get_node(tn, "mapped_inlet");
312 cs_gui_node_get_status_int(tn, &mapped_inlet);
313
314 if (mapped_inlet) {
315 cs_real_t coord_shift[3] = {0., 0., 0.};
316 const char *tname[] = {"translation_x",
317 "translation_y",
318 "translation_z"};
319
320 for (int i = 0; i < 3; i++) {
321 cs_tree_node_t *node = cs_tree_get_node(tn, tname[i]);
322
323 const cs_real_t *v = NULL;
324 v = cs_tree_node_get_values_real(node);
325 if (v != NULL )
326 coord_shift[i] = v[0];
327 }
328
329 cs_boundary_conditions_add_map(z->location_id,
330 CS_MESH_LOCATION_CELLS,
331 coord_shift,
332 0.1);
333 }
334 }
335
336 /*-----------------------------------------------------------------------------
337 * Value of velocity for sliding wall.
338 *
339 * parameters:
340 * tn_vp <-- tree node associated with BC velocity/pressure
341 * z_name <-- zone name
342 *----------------------------------------------------------------------------*/
343
344 static void
_sliding_wall(cs_tree_node_t * tn_vp,const char * z_name)345 _sliding_wall(cs_tree_node_t *tn_vp,
346 const char *z_name)
347 {
348 const char f_name[] = "velocity";
349 cs_field_t *f = cs_field_by_name(f_name);
350
351 cs_real_t value[3] = {0, 0, 0};
352
353 for (cs_tree_node_t *tn = cs_tree_node_get_child(tn_vp, "dirichlet");
354 tn != NULL;
355 tn = cs_tree_node_get_next_of_name(tn)) {
356 const char *name = cs_gui_node_get_tag(tn, "name");
357 int c_id = -1;
358 cs_gui_node_get_child_int(tn, "component", &c_id);
359 if (strcmp("velocity", name) == 0 && c_id > -1 && c_id < f->dim) {
360 const cs_real_t *v = cs_tree_node_get_values_real(tn);
361 if (v != NULL)
362 value[c_id] = v[0];
363 }
364 }
365
366 cs_equation_add_bc_by_value(_get_equation_param(f_name),
367 CS_PARAM_BC_DIRICHLET,
368 z_name,
369 value);
370 }
371
372 /*-----------------------------------------------------------------------------
373 * Values for turbulence variable for the current inlet.
374 *
375 * parameters:
376 * tn_bc <-- tree node associated with inlet BC
377 * izone <-- associated BC zone id
378 *----------------------------------------------------------------------------*/
379
380 static void
_inlet_turbulence(cs_tree_node_t * tn_bc,int izone)381 _inlet_turbulence(cs_tree_node_t *tn_bc,
382 int izone)
383 {
384 cs_tree_node_t *tn_t = cs_tree_node_get_child(tn_bc, "turbulence");
385 const char *choice = cs_tree_node_get_tag(tn_t, "choice");
386 if (choice == NULL)
387 return;
388
389 if (cs_gui_strcmp(choice, "hydraulic_diameter"))
390 boundaries->icalke[izone] = 1;
391 else if(cs_gui_strcmp(choice, "turbulent_intensity"))
392 boundaries->icalke[izone] = 2;
393 else if(cs_gui_strcmp(choice, "formula")) {
394 boundaries->icalke[izone] = 0;
395 return;
396 }
397 else
398 return;
399
400 cs_gui_node_get_child_real(tn_t, "hydraulic_diameter",
401 &boundaries->dh[izone]);
402
403 if (cs_gui_strcmp(choice, "turbulent_intensity")) {
404 const cs_real_t *v
405 = cs_tree_node_get_child_values_real(tn_t, "turbulent_intensity");
406 if (v != NULL)
407 boundaries->xintur[izone] = v[0] * 0.01;
408 }
409 }
410
411 /*-----------------------------------------------------------------------------
412 * get scalar's values
413 *
414 * parameters:
415 * tn_bc <-- tree node associated with boundary condition
416 * izone <-- associated zone id
417 * f_id <-- field id
418 *----------------------------------------------------------------------------*/
419
420 static void
_boundary_scalar(cs_tree_node_t * tn_bc,int izone,int f_id)421 _boundary_scalar(cs_tree_node_t *tn_bc,
422 int izone,
423 int f_id)
424 {
425 cs_field_t *f = cs_field_by_id(f_id);
426 const int dim = f->dim;
427
428 cs_tree_node_t *tn_s = cs_tree_node_get_child(tn_bc, "scalar");
429 tn_s = cs_tree_node_get_sibling_with_tag(tn_s, "name", f->name);
430
431 if (dim > 1)
432 tn_s = cs_tree_node_get_child(tn_s, "component");
433
434 cs_equation_param_t *eqp = _get_equation_param(f->name);
435 const char *z_name = boundaries->label[izone];
436 const char *choice = cs_tree_node_get_tag(tn_s, "choice");
437 const char *cnv = cs_tree_node_get_tag(tn_s, "convert");
438
439 if (cnv != NULL) {
440 if (f == CS_F_(h) && strcmp(cnv, "temperature") == 0)
441 boundaries->t_to_h[izone] = true;
442 else
443 bft_error
444 (__FILE__, __LINE__, 0,
445 _("%s: conversion for field %s from variable %s not handled."),
446 __func__, f->name, cnv);
447 }
448
449 cs_real_t value[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
450 assert(dim <= 9);
451
452 bool possibly_incomplete = false;
453
454 /* FIXME: we should not need a loop over components, but
455 directly use vector values; if we do not yet have
456 multidimensional user variables in the GUI, we can handle
457 this more cleanly */
458
459 for (int i = 0; i < dim; i++) {
460
461 /* All components should use the same BC type */
462
463 if (i > 0 && choice != NULL) {
464 const char *choice_c = cs_tree_node_get_tag(tn_s, "choice");
465 if (choice_c != NULL)
466 if (strcmp(choice, choice_c))
467 bft_error
468 (__FILE__, __LINE__, 0,
469 _("%s: for field %s on zone %s,\n"
470 "BC types are mismatched (%s on component 0, %s on component %d."),
471 __func__, f->name, z_name, choice, choice_c, i);
472 }
473
474 if (choice != NULL) {
475
476 if (! strcmp(choice, "dirichlet")) {
477 const cs_real_t *v = cs_tree_node_get_child_values_real(tn_s, choice);
478 if (v != NULL) {
479 value[i] = *v;
480 }
481 else
482 possibly_incomplete = true;
483
484 /* T to H conversion not handled using xdef yet. */
485 if (boundaries->t_to_h[izone]) {
486 boundaries->type_code[f_id][izone] = DIRICHLET_CNV;
487 boundaries->values[f_id][izone * dim + i].val1 = v[0];
488 }
489 }
490 else if (! strcmp(choice, "neumann")) {
491 const cs_real_t *v = cs_tree_node_get_child_values_real(tn_s, choice);
492 if (v != NULL) {
493 value[i] = *v;
494 }
495 }
496 else if (! strcmp(choice, "dirichlet_formula")) {
497
498 const char *s = cs_tree_node_get_child_value_str(tn_s, choice);
499 if (s != NULL) {
500 boundaries->type_code[f_id][izone] = DIRICHLET_FORMULA;
501 boundaries->scalar_e[f_id][izone * dim + i] = true;
502 }
503
504 }
505 else if (! strcmp(choice, "neumann_formula")) {
506
507 const char *s = cs_tree_node_get_child_value_str(tn_s, choice);
508 if (s != NULL) {
509 boundaries->type_code[f_id][izone] = NEUMANN_FORMULA;
510 boundaries->scalar_e[f_id][izone * dim + i] = true;
511 }
512 }
513 else if (! strcmp(choice, "exchange_coefficient_formula")) {
514
515 const char *s = cs_tree_node_get_child_value_str(tn_s, choice);
516 if (s != NULL) {
517 boundaries->type_code[f_id][izone] = EXCHANGE_COEFF_FORMULA;
518 boundaries->scalar_e[f_id][izone * dim + i] = true;
519 }
520 }
521 else if (! strcmp(choice, "exchange_coefficient")) {
522
523 const cs_real_t *v;
524
525 v = cs_tree_node_get_child_values_real(tn_s, "dirichlet");
526 if (v != NULL)
527 boundaries->values[f_id][izone * dim + i].val1 = v[0];
528
529 v = cs_tree_node_get_child_values_real(tn_s, "exchange_coefficient");
530 if (v != NULL) {
531 boundaries->type_code[f_id][izone] = EXCHANGE_COEFF;
532 boundaries->values[f_id][izone * dim + i].val2 = v[0];
533 }
534
535 }
536 else if (cs_gui_strcmp(choice, "dirichlet_implicit")) {
537 boundaries->type_code[f_id][izone] = DIRICHLET_IMPLICIT;
538 }
539 else if (cs_gui_strcmp(choice, "neumann_implicit")) {
540 boundaries->type_code[f_id][izone] = NEUMANN_IMPLICIT;
541 }
542 }
543
544 if (f->dim > 1)
545 tn_s = cs_tree_node_get_next_of_name(tn_s);
546 }
547
548 /* Now define appropriate equation parameters */
549
550 if (choice != NULL && cnv == NULL) {
551
552 if (! strcmp(choice, "dirichlet")) {
553
554 /* Some specific models may have set value already, so
555 if the value here is the default, it should probably be
556 ignored (the XML/tree structure should be improved to
557 avoid this type of problem )*/
558 if ( possibly_incomplete == false
559 || cs_equation_find_bc(eqp, z_name) == NULL)
560 cs_equation_add_bc_by_value(eqp,
561 CS_PARAM_BC_DIRICHLET,
562 z_name,
563 value);
564 }
565
566 else if (! strcmp(choice, "neumann"))
567 cs_equation_add_bc_by_value(eqp,
568 CS_PARAM_BC_NEUMANN,
569 z_name,
570 value);
571
572 }
573 }
574
575 /*-----------------------------------------------------------------------------
576 * Get coal's data for inlet. Check if the current zone is an inlet only
577 * for an oxydant, or for oxydant and coal.
578 *
579 * parameters:
580 * tn_vp <-- tree node associated with velocity and pressure
581 * izone <-- associated zone id
582 *----------------------------------------------------------------------------*/
583
584 static void
_inlet_coal(cs_tree_node_t * tn_vp,int izone)585 _inlet_coal(cs_tree_node_t *tn_vp,
586 int izone)
587 {
588 const int n_coals = cs_glob_combustion_model->coal.n_coals;
589 const int *nclpch = cs_glob_combustion_model->coal.n_classes_per_coal;
590
591 int _n_coals = 0;
592
593 /* Count coal definitions */
594
595 for (cs_tree_node_t *tn0 = cs_tree_node_get_child(tn_vp, "coal");
596 tn0 != NULL;
597 tn0 = cs_tree_node_get_next_of_name(tn0), _n_coals++) {
598
599 const char *name = cs_tree_node_get_tag(tn0, "name");
600 if (name == NULL)
601 continue;
602
603 int icoal = -1;
604 if (sscanf(name + strlen("coal"), "%d", &icoal) != 1)
605 continue;
606 icoal -= 1;
607 if (icoal +1 > n_coals)
608 continue;
609
610 /* mass flow rate of coal */
611 const cs_real_t *v = cs_tree_node_get_child_values_real(tn0, "flow1");
612 if (v != NULL)
613 boundaries->qimpcp[izone][icoal] = v[0];
614
615 /* temperature of coal */
616 v = cs_tree_node_get_child_values_real(tn0, "temperature");
617 if (v != NULL)
618 boundaries->timpcp[izone][icoal] = v[0];
619
620 /* loop on number of class by coal for ratio (%) stored in distch */
621
622 for (int iclass = 0; iclass < nclpch[icoal]; iclass++) {
623
624 char classname[32];
625 snprintf(classname, 31, "class%2.2i", iclass+1);
626
627 cs_tree_node_t *tn1
628 = cs_tree_get_node_with_tag(tn0, "ratio", "name", classname);
629 v = cs_tree_node_get_values_real(tn1);
630 if (v != NULL)
631 boundaries->distch[izone][icoal][iclass] = v[0];
632
633 }
634
635 }
636
637 /* if there is no coal, it is an inlet only for oxydant */
638 if (_n_coals == 0) {
639 boundaries->ientat[izone] = 1;
640 boundaries->ientcp[izone] = 0;
641 }
642 else {
643 boundaries->ientat[izone] = 0;
644 boundaries->ientcp[izone] = 1;
645 if (_n_coals != n_coals)
646 bft_error(__FILE__, __LINE__, 0,
647 _("Invalid number of coal: %i xml: %i\n"),
648 n_coals, _n_coals);
649 }
650 }
651
652 /*-----------------------------------------------------------------------------
653 * Get gas combustion's data for inlet.
654 *
655 * Check if the current zone is an inlet only for an oxydant,
656 * or for oxydant and coal.
657 *
658 * parameters:
659 * tn_vp <-- tree node associated with velocity and pressure
660 * izone <-- associated zone id
661 *----------------------------------------------------------------------------*/
662
663 static void
_inlet_gas(cs_tree_node_t * tn_vp,int izone)664 _inlet_gas(cs_tree_node_t *tn_vp,
665 int izone)
666 {
667 cs_tree_node_t *tn = cs_tree_get_node(tn_vp, "gas_type");
668 const char *choice = cs_gui_node_get_tag(tn, "choice");
669
670 if (cs_gui_strcmp(choice, "oxydant"))
671 boundaries->ientox[izone] = 1;
672
673 else if (cs_gui_strcmp(choice, "fuel"))
674 boundaries->ientfu[izone] = 1;
675
676 else if (cs_gui_strcmp(choice, "unburned")) {
677
678 boundaries->ientgf[izone] = 1;
679
680 cs_gui_node_get_child_real(tn_vp, "temperature",
681 &boundaries->tkent[izone]);
682 cs_gui_node_get_child_real(tn_vp, "fraction",
683 &boundaries->fment[izone]);
684
685 }
686 else if (cs_gui_strcmp(choice, "burned")) {
687
688 boundaries->ientgb[izone] = 1;
689
690 cs_gui_node_get_child_real(tn_vp, "temperature",
691 &boundaries->tkent[izone]);
692 cs_gui_node_get_child_real(tn_vp, "fraction",
693 &boundaries->fment[izone]);
694
695 }
696 }
697
698 /*-----------------------------------------------------------------------------
699 * Get compressible data for inlet.
700 *
701 * parameters:
702 * tn_vp <-- tree node associated with velocity and pressure
703 * izone <-- associated zone id
704 *----------------------------------------------------------------------------*/
705
706 static void
_inlet_compressible(cs_tree_node_t * tn_vp,int izone)707 _inlet_compressible(cs_tree_node_t *tn_vp,
708 int izone)
709 {
710 const cs_zone_t *z = cs_boundary_zone_by_id(izone + 1);
711
712 bool status;
713
714 cs_tree_node_t *tn = cs_tree_get_node(tn_vp, "compressible_type");
715 const char *choice = cs_gui_node_get_tag(tn, "choice");
716
717 if (cs_gui_strcmp(choice, "imposed_inlet")) {
718
719 cs_real_t te_in = cs_math_infinite_r;
720
721 boundaries->itype[izone] = CS_ESICF;
722
723 tn = cs_tree_node_get_child(tn_vp, "pressure");
724 status = false;
725 cs_gui_node_get_status_bool(tn, &status);
726 if (status)
727 cs_gui_node_get_real(tn, &boundaries->prein[izone]);
728
729 tn = cs_tree_node_get_child(tn_vp, "density");
730 status = false;
731 cs_gui_node_get_status_bool(tn, &status);
732 if (status)
733 cs_gui_node_get_real(tn, &boundaries->rhoin[izone]);
734
735 tn = cs_tree_node_get_child(tn_vp, "temperature");
736 status = false;
737 cs_gui_node_get_status_bool(tn, &status);
738 if (status)
739 cs_gui_node_get_real(tn, &boundaries->tempin[izone]);
740
741 tn = cs_tree_node_get_child(tn_vp, "energy");
742 status = false;
743 cs_gui_node_get_status_bool(tn, &status);
744 if (status)
745 cs_gui_node_get_real(tn, &te_in);
746
747 cs_equation_param_t *eqp = _get_equation_param("total_energy");
748 cs_equation_remove_bc(eqp, z->name);
749 cs_equation_add_bc_by_value(eqp,
750 CS_PARAM_BC_DIRICHLET,
751 z->name,
752 &te_in);
753
754 }
755 else if (cs_gui_strcmp(choice, "subsonic_inlet_PH")) {
756
757 boundaries->itype[izone] = CS_EPHCF;
758
759 cs_gui_node_get_child_real
760 (tn_vp, "total_pressure", &boundaries->prein[izone]);
761
762 cs_real_t h_in = cs_math_infinite_r;
763 cs_gui_node_get_child_real(tn_vp, "enthalpy", &h_in);
764
765 cs_equation_param_t *eqp = _get_equation_param("total_energy");
766 cs_equation_remove_bc(eqp, z->name);
767 cs_equation_add_bc_by_value(eqp,
768 CS_PARAM_BC_DIRICHLET,
769 z->name,
770 &h_in);
771
772 }
773 }
774
775 /*-----------------------------------------------------------------------------
776 * Get compressible data for inlet.
777 *
778 * parameters:
779 * tn_bc <-- tree node associated with boundary zone
780 * izone <-- associated zone id
781 *----------------------------------------------------------------------------*/
782
783 static void
_outlet_compressible(cs_tree_node_t * tn_bc,int izone)784 _outlet_compressible(cs_tree_node_t *tn_bc,
785 int izone)
786 {
787 const char *z_name = boundaries->label[izone];
788
789 cs_tree_node_t *tn = cs_tree_node_get_child(tn_bc, "compressible_type");
790
791 const char *choice = cs_tree_node_get_tag(tn, "choice");
792
793 if (cs_gui_strcmp(choice, "supersonic_outlet")) {
794 boundaries->itype[izone] = CS_SSPCF;
795 }
796 else if (cs_gui_strcmp(choice, "subsonic_outlet")) {
797 boundaries->itype[izone] = CS_SOPCF;
798
799 const char name[] = "pressure";
800 tn = cs_tree_node_get_child(tn_bc, "dirichlet");
801 tn = cs_tree_node_get_sibling_with_tag(tn, "name", name);
802
803 if (tn != NULL) {
804 cs_real_t value = 0;
805 const cs_real_t *v = cs_tree_node_get_values_real(tn);
806 if (v != NULL)
807 value = v[0];
808 cs_equation_add_bc_by_value(_get_equation_param(name),
809 CS_PARAM_BC_DIRICHLET,
810 z_name,
811 &value);
812 }
813 }
814 }
815
816 /*-----------------------------------------------------------------------------
817 * Get pressure value for darcy (inlet/outlet/groundwater).
818 *
819 * parameters:
820 * tn_bc <-- tree node associated with boundary conditions
821 * izone <-- associated zone id
822 *----------------------------------------------------------------------------*/
823
824 static void
_boundary_darcy(cs_tree_node_t * tn_bc,int izone)825 _boundary_darcy(cs_tree_node_t *tn_bc,
826 int izone)
827 {
828 const char *z_name = boundaries->label[izone];
829
830 cs_tree_node_t *tn = cs_tree_node_get_child(tn_bc, "hydraulicHead");
831 const char *choice = cs_gui_node_get_tag(tn, "choice");
832
833 tn = cs_tree_node_get_child(tn_bc, choice);
834 tn = cs_tree_node_get_sibling_with_tag(tn, "name", "hydraulic_head");
835
836 cs_equation_param_t *eqp = cs_field_get_equation_param(CS_F_(head));
837 if (eqp == NULL)
838 eqp = _get_equation_param("pressure_head"); /* CDO version */
839
840 if (cs_gui_strcmp(choice, "dirichlet")) {
841 cs_real_t value = 0;
842 cs_gui_node_get_real(tn, &value);
843 cs_equation_add_bc_by_value(eqp,
844 CS_PARAM_BC_DIRICHLET,
845 z_name,
846 &value);
847 }
848 else if (cs_gui_strcmp(choice, "neumann")) {
849 /* Vector values per component for CDO, scalar (1st component) for legacy */
850 cs_real_t value[3] = {0, 0, 0};
851 cs_gui_node_get_real(tn, value);
852 cs_equation_add_bc_by_value(eqp,
853 CS_PARAM_BC_NEUMANN,
854 z_name,
855 value);
856 }
857 else if (cs_gui_strcmp(choice, "dirichlet_formula")) {
858 if (tn == NULL) { /* compatibility with inconsistant tag */
859 tn = cs_tree_node_get_child(tn_bc, choice);
860 tn = cs_tree_node_get_sibling_with_tag(tn, "name", "hydraulicHead");
861 }
862 const char *formula = cs_tree_node_get_child_value_str(tn, "formula");
863 if (formula != NULL)
864 boundaries->groundwat_e[izone] = true;
865 else {
866 bft_printf("Warning : groundwater flow boundary conditions\n"
867 " without formula for hydraulic head.\n");
868 }
869 }
870 }
871
872 /*-----------------------------------------------------------------------------
873 * Get pressure value for imposed pressure boundary
874 *
875 * parameters:
876 * tn_bc <-- tree node associated with boundary conditions
877 * z_name <-- id of the current zone
878 *----------------------------------------------------------------------------*/
879
880 static void
_boundary_imposed_pressure(cs_tree_node_t * tn_bc,const char * z_name)881 _boundary_imposed_pressure(cs_tree_node_t *tn_bc,
882 const char *z_name)
883 {
884 const char name[] = "pressure";
885 cs_tree_node_t *tn = cs_tree_node_get_child(tn_bc, "dirichlet");
886 tn = cs_tree_node_get_sibling_with_tag(tn, "name", name);
887
888 cs_real_t value = 0;
889 cs_gui_node_get_real(tn, &value);
890
891 cs_equation_param_t *eqp = _get_equation_param(name);
892 cs_equation_add_bc_by_value(eqp,
893 CS_PARAM_BC_DIRICHLET,
894 z_name,
895 &value);
896 }
897
898 /*-----------------------------------------------------------------------------
899 * Helper to get the face list for the izone
900 *
901 * parameters:
902 * label <-- boundary label
903 * n_faces --> number of faces
904 *
905 * returns:
906 * pointer to face list
907 *----------------------------------------------------------------------------*/
908
909 static const cs_lnum_t *
_get_boundary_faces(const char * label,cs_lnum_t * n_faces)910 _get_boundary_faces(const char *label,
911 cs_lnum_t *n_faces)
912 {
913 const cs_lnum_t *face_ids = NULL;
914
915 const cs_zone_t *z = cs_boundary_zone_by_name(label);
916
917 *n_faces = z->n_elts;
918 face_ids = z->elt_ids;
919
920 return face_ids;
921 }
922
923 /*----------------------------------------------------------------------------
924 * Boundary conditions treatment: global structure initialization
925 *----------------------------------------------------------------------------*/
926
927 static void
_init_boundaries(void)928 _init_boundaries(void)
929 {
930 int icharb, iclass;
931
932 assert(boundaries == NULL);
933 int n_fields = cs_field_n_fields();
934
935 int n_zones = cs_tree_get_node_count(cs_glob_tree,
936 "boundary_conditions/boundary");
937
938 bool solid_fuels = false;
939 if ( cs_glob_physical_model_flag[CS_COMBUSTION_PCLC] > -1
940 || cs_glob_physical_model_flag[CS_COMBUSTION_COAL] > -1)
941 solid_fuels = true;
942 bool gas_combustion = false;
943 for (cs_physical_model_type_t m_type = CS_COMBUSTION_3PT;
944 m_type <= CS_COMBUSTION_LW;
945 m_type++) {
946 if (cs_glob_physical_model_flag[m_type] > -1)
947 gas_combustion = true;
948 }
949 if ( cs_glob_physical_model_flag[CS_COMBUSTION_PCLC] > -1
950 || cs_glob_physical_model_flag[CS_COMBUSTION_COAL] > -1)
951 solid_fuels = true;
952
953 BFT_MALLOC(boundaries, 1, cs_gui_boundary_t);
954
955 boundaries->n_fields = n_fields;
956 boundaries->n_zones = n_zones;
957
958 BFT_MALLOC(boundaries->label, n_zones, const char *);
959 BFT_MALLOC(boundaries->nature, n_zones, const char *);
960 BFT_MALLOC(boundaries->bc_num, n_zones, int);
961
962 BFT_MALLOC(boundaries->iqimp, n_zones, int);
963
964 boundaries->ientfu = NULL;
965 boundaries->ientox = NULL;
966 boundaries->ientgb = NULL;
967 boundaries->ientgf = NULL;
968 boundaries->ientat = NULL;
969 boundaries->ientcp = NULL;
970
971 BFT_MALLOC(boundaries->icalke, n_zones, int);
972 BFT_MALLOC(boundaries->qimp, n_zones, double);
973
974 boundaries->inmoxy = NULL;
975 boundaries->timpat = NULL;
976 boundaries->tkent = NULL;
977 boundaries->qimpcp = NULL;
978 boundaries->timpcp = NULL;
979 boundaries->fment = NULL;
980 boundaries->itype = NULL;
981 boundaries->prein = NULL;
982 boundaries->rhoin = NULL;
983 boundaries->tempin = NULL;
984
985 BFT_MALLOC(boundaries->dh, n_zones, double);
986 BFT_MALLOC(boundaries->xintur, n_zones, double);
987 BFT_MALLOC(boundaries->type_code, n_fields, int *);
988 BFT_MALLOC(boundaries->values, n_fields, cs_val_t *);
989
990 boundaries->distch = NULL;
991
992 BFT_MALLOC(boundaries->rough, n_zones, double);
993 BFT_MALLOC(boundaries->norm, n_zones, double);
994 BFT_MALLOC(boundaries->dir, n_zones, cs_real_3_t);
995
996 BFT_MALLOC(boundaries->t_to_h, n_zones, bool);
997 BFT_MALLOC(boundaries->velocity_e, n_zones, bool);
998 BFT_MALLOC(boundaries->direction_e, n_zones, bool);
999 BFT_MALLOC(boundaries->scalar_e, n_fields, bool *);
1000 BFT_MALLOC(boundaries->head_loss_e, n_zones, bool);
1001
1002 boundaries->groundwat_e = NULL;
1003 boundaries->meteo = NULL;
1004
1005 if (solid_fuels) {
1006
1007 const cs_combustion_model_t *cm = cs_glob_combustion_model;
1008
1009 BFT_MALLOC(boundaries->ientat, n_zones, int);
1010 BFT_MALLOC(boundaries->inmoxy, n_zones, int);
1011 BFT_MALLOC(boundaries->timpat, n_zones, double);
1012 BFT_MALLOC(boundaries->ientcp, n_zones, int);
1013 BFT_MALLOC(boundaries->qimpcp, n_zones, double *);
1014 BFT_MALLOC(boundaries->timpcp, n_zones, double *);
1015 BFT_MALLOC(boundaries->distch, n_zones, double **);
1016
1017 for (int izone=0; izone < n_zones; izone++) {
1018 BFT_MALLOC(boundaries->qimpcp[izone], cm->coal.n_coals, double);
1019 BFT_MALLOC(boundaries->timpcp[izone], cm->coal.n_coals, double);
1020 BFT_MALLOC(boundaries->distch[izone], cm->coal.n_coals, double *);
1021
1022 for (icharb = 0; icharb < cm->coal.n_coals; icharb++)
1023 BFT_MALLOC(boundaries->distch[izone][icharb],
1024 cm->coal.n_classes_per_coal[icharb],
1025 double);
1026 }
1027 }
1028 else if (gas_combustion) {
1029 BFT_MALLOC(boundaries->ientfu, n_zones, int);
1030 BFT_MALLOC(boundaries->ientox, n_zones, int);
1031 BFT_MALLOC(boundaries->ientgb, n_zones, int);
1032 BFT_MALLOC(boundaries->ientgf, n_zones, int);
1033 BFT_MALLOC(boundaries->tkent, n_zones, double);
1034 BFT_MALLOC(boundaries->fment, n_zones, double);
1035 }
1036 else if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
1037 BFT_MALLOC(boundaries->itype, n_zones, int);
1038 BFT_MALLOC(boundaries->prein, n_zones, double);
1039 BFT_MALLOC(boundaries->rhoin, n_zones, double);
1040 BFT_MALLOC(boundaries->tempin, n_zones, double);
1041 }
1042 else if (cs_glob_physical_model_flag[CS_GROUNDWATER] > -1) {
1043 BFT_MALLOC(boundaries->groundwat_e, n_zones, bool);
1044 }
1045
1046 if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] > -1)
1047 BFT_MALLOC(boundaries->meteo, n_zones, cs_meteo_t);
1048 else
1049 boundaries->meteo = NULL;
1050
1051 for (int f_id = 0; f_id < n_fields; f_id++) {
1052 const cs_field_t *f = cs_field_by_id(f_id);
1053
1054 if (f->type & CS_FIELD_VARIABLE) {
1055 BFT_MALLOC(boundaries->type_code[f->id], n_zones, int);
1056 BFT_MALLOC(boundaries->values[f->id], n_zones * f->dim, cs_val_t);
1057 BFT_MALLOC(boundaries->scalar_e[f->id], n_zones * f->dim, bool);
1058 }
1059 }
1060
1061 /* initialize for each zone */
1062 for (int izone = 0; izone < n_zones; izone++) {
1063 boundaries->iqimp[izone] = 0;
1064 boundaries->qimp[izone] = 0;
1065 boundaries->norm[izone] = 0;
1066 boundaries->icalke[izone] = 0;
1067 boundaries->dh[izone] = 0;
1068 boundaries->xintur[izone] = 0;
1069 boundaries->rough[izone] = -999;
1070 boundaries->t_to_h[izone] = false;
1071 boundaries->velocity_e[izone] = false;
1072 boundaries->direction_e[izone] = false;
1073 boundaries->head_loss_e[izone] = false;
1074
1075 if (solid_fuels) {
1076 const cs_combustion_model_t *cm = cs_glob_combustion_model;
1077
1078 boundaries->ientat[izone] = 0;
1079 boundaries->inmoxy[izone] = 1;
1080 boundaries->ientcp[izone] = 0;
1081 boundaries->timpat[izone] = 0;
1082
1083 for (icharb = 0; icharb < cm->coal.n_coals; icharb++) {
1084 boundaries->qimpcp[izone][icharb] = 0;
1085 boundaries->timpcp[izone][icharb] = 0;
1086
1087 for (iclass = 0; iclass < cm->coal.n_classes_per_coal[icharb]; iclass++)
1088 boundaries->distch[izone][icharb][iclass] = 0;
1089 }
1090 }
1091
1092 else if (gas_combustion) {
1093 boundaries->ientfu[izone] = 0;
1094 boundaries->ientox[izone] = 0;
1095 boundaries->ientgb[izone] = 0;
1096 boundaries->ientgf[izone] = 0;
1097 boundaries->tkent[izone] = 0;
1098 boundaries->fment[izone] = 0;
1099 }
1100
1101 else if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
1102 boundaries->itype[izone] = 0;
1103 boundaries->prein[izone] = cs_math_infinite_r;
1104 boundaries->rhoin[izone] = 0;
1105 boundaries->tempin[izone] = cs_math_infinite_r;
1106 }
1107
1108 else if (cs_glob_physical_model_flag[CS_GROUNDWATER] > -1) {
1109 boundaries->groundwat_e[izone] = false;
1110 }
1111
1112 else if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] > -1) {
1113 boundaries->meteo[izone].read_data = 0;
1114 boundaries->meteo[izone].automatic = 0;
1115 }
1116 }
1117
1118 /* Initialization of boundary->type_code and boundary->values */
1119 for (int f_id = 0; f_id < n_fields; f_id++) {
1120 const cs_field_t *f = cs_field_by_id(f_id);
1121
1122 if (f->type & CS_FIELD_VARIABLE) {
1123 int i = f->id;
1124 for (int izone = 0; izone < n_zones; izone++) {
1125 boundaries->type_code[i][izone] = -1;
1126 for (int ii = 0; ii < f->dim; ii++) {
1127 boundaries->values[i][izone * f->dim + ii].val1 = 1.e30;
1128 boundaries->values[i][izone * f->dim + ii].val2 = 1.e30;
1129 boundaries->scalar_e[i][izone * f->dim + ii] = false;
1130 }
1131 }
1132 }
1133 }
1134
1135 /* filling of the "boundaries" structure */
1136
1137 cs_tree_node_t *tn_b0 = cs_tree_get_node(cs_glob_tree,
1138 "boundary_conditions");
1139
1140 int izone = 0;
1141
1142 /* Complete boundary zone definitions */
1143
1144 for (cs_tree_node_t *tn = cs_tree_get_node(tn_b0, "boundary");
1145 tn != NULL;
1146 tn = cs_tree_node_get_next_of_name(tn), izone++) {
1147
1148 /* nature, label and description of the ith boundary zone;
1149 zones are shifted by 1, as default zone 0 is defined
1150 first, before GUI-based definitions (and non-GUI-based
1151 user definitions come last). */
1152
1153 const char *label = cs_tree_node_get_tag(tn, "label");
1154 const char *nature = cs_tree_node_get_tag(tn, "nature");
1155
1156 int bc_num = izone+1;
1157
1158 const int *vi = cs_tree_node_get_child_values_int(tn, "name");
1159 if (vi != NULL)
1160 bc_num = vi[0];
1161
1162 /* label of the ith initialization zone */
1163
1164 const cs_zone_t *z = cs_boundary_zone_by_id(izone + 1);
1165
1166 assert(strcmp(label, z->name) == 0);
1167
1168 /* Note: boundary nature is determined again later for most cases,
1169 but at least symmetry is only defined here, and ALE can define
1170 "wall" sections instead of the appropriate type to handle
1171 fixed sections, so also predefine the boundary nature here */
1172
1173 boundaries->label[izone] = z->name;
1174 boundaries->nature[izone] = nature;
1175 boundaries->bc_num[izone] = bc_num;
1176
1177 }
1178
1179 cs_wall_functions_t *wall_fnt = cs_get_glob_wall_functions();
1180
1181 /* Now loop on boundary condition definitions proper */
1182
1183 cs_tree_node_t *tn_b1 = (tn_b0 != NULL) ? tn_b0->children : tn_b0;
1184
1185 for (cs_tree_node_t *tn = tn_b1; tn != NULL; tn = tn->next) {
1186
1187 if (cs_gui_strcmp(tn->name, "boundary")) /* handled in previous loop */
1188 continue;
1189
1190 const cs_zone_t *z = NULL;
1191
1192 const char *nature = tn->name;
1193 const char *label = cs_tree_node_get_tag(tn, "label");
1194 if (label != NULL)
1195 z = cs_boundary_zone_by_name_try(label);
1196
1197 if (z == NULL) /* may occur when "dead" leaves or present in tree */
1198 continue;
1199
1200 izone = z->id - 1;
1201 assert(izone >= 0);
1202
1203 /* Note: ALE may define mesh BCs as "wall" zones even where this is not
1204 appropriate, so skip it here (as this is handled elsewhere) */
1205 if (strcmp(boundaries->nature[izone], nature))
1206 continue;
1207
1208 if (cs_gui_strcmp(nature, "inlet")) {
1209
1210 _check_and_add_mapped_inlet(label, z);
1211
1212 cs_tree_node_t *tn_vp
1213 = cs_tree_node_get_child(tn, "velocity_pressure");
1214
1215 if (cs_glob_physical_model_flag[CS_GROUNDWATER] < 0) {
1216
1217 const char *choice_v = cs_gui_node_get_tag(tn_vp, "choice");
1218 const char *choice_d = cs_gui_node_get_tag(tn_vp, "direction");
1219
1220 /* Inlet: velocity */
1221
1222 if (cs_gui_strcmp(choice_v, "norm")) {
1223 cs_gui_node_get_child_real
1224 (tn_vp, "norm", &boundaries->norm[izone]);
1225 }
1226 else if (cs_gui_strcmp(choice_v, "flow1")) {
1227 cs_gui_node_get_child_real
1228 (tn_vp, "flow1", &boundaries->qimp[izone]);
1229 boundaries->iqimp[izone] = 1;
1230 }
1231 else if (cs_gui_strcmp(choice_v, "flow2")) {
1232 cs_gui_node_get_child_real
1233 (tn_vp, "flow2", &boundaries->qimp[izone]);
1234 boundaries->iqimp[izone] = 2;
1235 }
1236 else if (cs_gui_strcmp(choice_v, "norm_formula")) {
1237 if (cs_tree_node_get_child_value_str(tn_vp, choice_v) != NULL)
1238 boundaries->velocity_e[izone] = true;
1239 }
1240 else if (cs_gui_strcmp(choice_v, "flow1_formula")) {
1241 if (cs_tree_node_get_child_value_str(tn_vp, choice_v) != NULL)
1242 boundaries->velocity_e[izone] = true;
1243 boundaries->iqimp[izone] = 1;
1244 }
1245 else if (cs_gui_strcmp(choice_v, "flow2_formula")) {
1246 if (cs_tree_node_get_child_value_str(tn_vp, choice_v) != NULL)
1247 boundaries->velocity_e[izone] = true;
1248 boundaries->iqimp[izone] = 2;
1249 }
1250
1251 if ( cs_gui_strcmp(choice_d, "coordinates")
1252 || cs_gui_strcmp(choice_d, "translation")) {
1253 cs_real_t *dir = boundaries->dir[izone];
1254 cs_gui_node_get_child_real(tn_vp, "direction_x", dir);
1255 cs_gui_node_get_child_real(tn_vp, "direction_y", dir+1);
1256 cs_gui_node_get_child_real(tn_vp, "direction_z", dir+2);
1257 }
1258 else if (cs_gui_strcmp(choice_d, "formula")) {
1259 if ( cs_tree_node_get_child_value_str(tn_vp, "direction_formula")
1260 != NULL)
1261 boundaries->direction_e[izone] = true;
1262 }
1263 }
1264
1265 /* Inlet: data for COAL COMBUSTION */
1266 if (solid_fuels) {
1267 cs_gui_node_get_child_real
1268 (tn_vp, "temperature", &boundaries->timpat[izone]);
1269 cs_gui_node_get_child_int
1270 (tn_vp, "oxydant", &boundaries->inmoxy[izone]);
1271 _inlet_coal(tn_vp, izone);
1272 }
1273
1274 /* Inlet: data for GAS COMBUSTION */
1275 if (gas_combustion)
1276 _inlet_gas(tn_vp, izone);
1277
1278 /* Inlet: data for COMPRESSIBLE MODEL */
1279 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1)
1280 _inlet_compressible(tn_vp, izone);
1281
1282 /* Inlet: data for ATMOSPHERIC FLOWS */
1283 if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] > -1) {
1284 if (cs_glob_atmo_option->meteo_profile > 0) {
1285 cs_gui_node_get_child_status_int
1286 (tn_vp, "meteo_data", &boundaries->meteo[izone].read_data);
1287 cs_gui_node_get_child_status_int
1288 (tn_vp, "meteo_automatic", &boundaries->meteo[izone].automatic);
1289 }
1290 }
1291
1292 /* Inlet: data for Darcy */
1293 if (cs_glob_physical_model_flag[CS_GROUNDWATER] > -1)
1294 _boundary_darcy(tn, izone);
1295
1296 /* Inlet: turbulence */
1297 _inlet_turbulence(tn, izone);
1298
1299 }
1300 else if (cs_gui_strcmp(nature, "wall")) {
1301
1302 /* sliding wall: velocity */
1303
1304 cs_tree_node_t *tn_vp
1305 = cs_tree_node_get_child(tn, "velocity_pressure");
1306
1307 if (tn_vp != NULL) {
1308 /* Wall: ROUGH */
1309 if ( wall_fnt->iwallf != CS_WALL_F_DISABLED
1310 && wall_fnt->iwallf != CS_WALL_F_1SCALE_POWER
1311 && wall_fnt->iwallf != CS_WALL_F_SCALABLE_2SCALES_LOG
1312 && wall_fnt->iwallf != CS_WALL_F_2SCALES_CONTINUOUS) {
1313 cs_gui_node_get_child_real(tn_vp, "roughness",
1314 &boundaries->rough[izone]);
1315 }
1316 }
1317 }
1318
1319 else if (cs_gui_strcmp(nature, "outlet")) {
1320 /* Outlet: data for ATMOSPHERIC FLOWS */
1321 if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] > -1) {
1322 _boundary_status_vp("outlet", label, "meteo_data",
1323 &boundaries->meteo[izone].read_data);
1324 _boundary_status_vp("outlet", label, "meteo_automatic",
1325 &boundaries->meteo[izone].automatic);
1326 }
1327
1328 /* Outlet: data for COMPRESSIBLE MODEL */
1329 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1)
1330 _outlet_compressible(tn, izone);
1331
1332 /* Inlet: data for darcy */
1333 if (cs_glob_physical_model_flag[CS_GROUNDWATER] > -1)
1334 _boundary_darcy(tn, izone);
1335 }
1336
1337 else if (cs_gui_strcmp(nature, "free_inlet_outlet")) {
1338 cs_tree_node_t *tn_hlf = cs_tree_get_node(tn, "headLoss/formula");
1339 const char *hl_formula = cs_tree_node_get_value_str(tn_hlf);
1340 if (hl_formula != NULL)
1341 boundaries->head_loss_e[izone] = true;
1342 else {
1343 bft_printf("Warning : free inlet outlet boundary conditions\n"
1344 " without external head loss definition\n");
1345 }
1346 }
1347 else if (cs_gui_strcmp(nature, "imposed_p_outlet")) {
1348 _boundary_imposed_pressure(tn, label);
1349 }
1350 else if (cs_gui_strcmp(nature, "groundwater")) {
1351 _boundary_darcy(tn, izone);
1352 }
1353
1354 /* for each zone */
1355 if (!cs_gui_strcmp(nature, "symmetry")) {
1356
1357 /* Thermal scalar */
1358
1359 cs_field_t *f_tm = cs_thermal_model_field();
1360
1361 if (f_tm != NULL) {
1362 if (boundaries->meteo == NULL)
1363 _boundary_scalar(tn, izone, f_tm->id);
1364 else if (boundaries->meteo[izone].read_data == 0)
1365 _boundary_scalar(tn, izone, f_tm->id);
1366 }
1367
1368 const char *scalar_sections[]
1369 = {"thermophysical_models/atmospheric_flows/variable",
1370 "thermophysical_models/joule_effect/variable"};
1371
1372 /* Meteo scalars only if required */
1373 if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] < 0)
1374 scalar_sections[0] = NULL;
1375 else {
1376 if (boundaries->meteo[izone].read_data != 0)
1377 scalar_sections[0] = NULL;
1378 }
1379
1380 /* Electric arc scalars only if required */
1381 if (cs_glob_physical_model_flag[CS_ELECTRIC_ARCS] < 0)
1382 scalar_sections[1] = NULL;
1383
1384 /* Loop on possible specific model scalar sections */
1385 for (int s_id = 0; s_id < 2; s_id++) {
1386 if (scalar_sections[s_id] == NULL)
1387 continue;
1388 for (cs_tree_node_t *tn_sv
1389 = cs_tree_get_node(cs_glob_tree, scalar_sections[s_id]);
1390 tn_sv != NULL;
1391 tn_sv = cs_tree_node_get_next_of_name(tn_sv)) {
1392 const char *_name = cs_gui_node_get_tag(tn_sv, "name");
1393 cs_field_t *f = cs_field_by_name_try(_name);
1394 if (f != NULL)
1395 _boundary_scalar(tn, izone, f->id);
1396 }
1397 }
1398
1399 /* User scalars */
1400 for (int f_id = 0; f_id < n_fields; f_id++) {
1401 const cs_field_t *f = cs_field_by_id(f_id);
1402 if ( (f->type & CS_FIELD_VARIABLE)
1403 && (f->type & CS_FIELD_USER))
1404 _boundary_scalar(tn, izone, f->id);
1405 }
1406 }
1407
1408 } /* for izones */
1409 }
1410
1411 /*----------------------------------------------------------------------------
1412 * Boundary conditions treatment: initialize and check zone info
1413 *
1414 * parameters:
1415 * n_b_faces <-- number of boundary faces
1416 * nozppm <-- max number of boundary conditions zone
1417 * izfppp <-- zone number for each boundary face
1418 *----------------------------------------------------------------------------*/
1419
1420 static void
_init_zones(const cs_lnum_t n_b_faces,const int * nozppm,int * izfppp)1421 _init_zones(const cs_lnum_t n_b_faces,
1422 const int *nozppm,
1423 int *izfppp)
1424 {
1425
1426 assert(boundaries != NULL);
1427
1428 int n_zones = cs_tree_get_node_count(cs_glob_tree,
1429 "boundary_conditions/boundary");
1430
1431 for (cs_lnum_t ifac = 0; ifac < n_b_faces; ifac++)
1432 izfppp[ifac] = 0;
1433
1434 for (int izone = 0; izone < n_zones; izone++) {
1435
1436 int zone_nbr = boundaries->bc_num[izone];
1437
1438 if (nozppm && zone_nbr > *nozppm)
1439 bft_error(__FILE__, __LINE__, 0,
1440 _("zone's label number %i is greater than %i,"
1441 " the maximum allowed \n"), zone_nbr, *nozppm);
1442
1443 cs_lnum_t n_faces = 0;
1444 const cs_lnum_t *face_ids
1445 = _get_boundary_faces(boundaries->label[izone], &n_faces);
1446
1447 /* check if faces are already marked with a zone number */
1448
1449 for (cs_lnum_t f_id = 0; f_id < n_faces; f_id++) {
1450 cs_lnum_t ifbr = face_ids[f_id];
1451 izfppp[ifbr] = zone_nbr;
1452 }
1453
1454 } /* for izone */
1455
1456 }
1457
1458 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1459
1460 /*============================================================================
1461 * Public Fortran function definitions
1462 *============================================================================*/
1463
1464 /*----------------------------------------------------------------------------
1465 * Boundary conditions treatment
1466 *
1467 * Remember: rdoccl[k][j][i] = rcodcl[ k * dim1 *dim2 + j *dim1 + i]
1468 *
1469 * Fortran Interface:
1470 *
1471 * subroutine uiclim
1472 * *****************
1473 *
1474 * integer nozppm <-- max number of boundary conditions zone
1475 * integer nclpch <-- number of simulated class per coals
1476 * integer iqimp <-- 1 if flow rate is applied
1477 * integer icalke <-- 1 for automatic turbulent boundary conditions
1478 * integer ientat <-- 1 for air temperature boundary conditions (coal)
1479 * integer ientcp <-- 1 for coal temperature boundary conditions (coal)
1480 * INTEGER INMOXY <-- coal: number of oxydant for the current inlet
1481 * integer ientox <-- 1 for an air fow inlet (gas combustion)
1482 * integer ientfu <-- 1 for fuel flow inlet (gas combustion)
1483 * integer ientgb <-- 1 for burned gas inlet (gas combustion)
1484 * integer ientgf <-- 1 for unburned gas inlet (gas combustion)
1485 * integer iprofm <-- atmospheric flows: on/off for profile from data
1486 * integer iautom <-- atmospheric flows: auto inlet/outlet flag
1487 * integer itypfb <-- type of boundary for each face
1488 * integer izfppp <-- zone number for each boundary face
1489 * integer icodcl <-- boundary conditions array type
1490 * double precision cgdfbo <-- boundary faces center of gravity
1491 * double precision qimp <-- inlet flow rate
1492 * double precision qimpat <-- inlet air flow rate (coal)
1493 * double precision qimpcp <-- inlet coal flow rate (coal)
1494 * double precision dh <-- hydraulic diameter
1495 * double precision xintur <-- turbulent intensity
1496 * double precision timpat <-- air temperature boundary conditions (coal)
1497 * double precision timpcp <-- inlet coal temperature (coal)
1498 * double precision tkent <-- inlet temperature (gas combustion)
1499 * double precision fment <-- Mean Mixture Fraction at Inlet (gas combustion)
1500 * double precision distch <-- ratio for each coal
1501 * integer nvar <-- dimension for rcodcl
1502 * double precision rcodcl <-- boundary conditions array value
1503 *----------------------------------------------------------------------------*/
1504
CS_PROCF(uiclim,UICLIM)1505 void CS_PROCF (uiclim, UICLIM)(const int *nozppm,
1506 int *iqimp,
1507 int *icalke,
1508 int *ientat,
1509 int *ientcp,
1510 int *inmoxy,
1511 int *ientox,
1512 int *ientfu,
1513 int *ientgf,
1514 int *ientgb,
1515 int *iprofm,
1516 int *iautom,
1517 int *itypfb,
1518 int *izfppp,
1519 int *icodcl,
1520 double *qimp,
1521 double *qimpat,
1522 double *qimpcp,
1523 double *dh,
1524 double *xintur,
1525 double *timpat,
1526 double *timpcp,
1527 double *tkent,
1528 double *fment,
1529 double *distch,
1530 int *nvar,
1531 double *rcodcl)
1532 {
1533 double norm = 0.;
1534
1535 const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
1536 const cs_lnum_t *b_face_cells = cs_glob_mesh->b_face_cells;
1537
1538 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
1539 const cs_real_t *b_face_surf = mq->b_face_surf;
1540 const cs_real_3_t *b_face_normal = (const cs_real_3_t *)mq->b_face_normal;
1541 const int n_fields = cs_field_n_fields();
1542
1543 const int ncharm = CS_COMBUSTION_MAX_COALS;
1544
1545 bool solid_fuels = false;
1546 if ( cs_glob_physical_model_flag[CS_COMBUSTION_PCLC] > -1
1547 || cs_glob_physical_model_flag[CS_COMBUSTION_COAL] > -1)
1548 solid_fuels = true;
1549 bool gas_combustion = false;
1550 for (cs_physical_model_type_t m_type = CS_COMBUSTION_3PT;
1551 m_type <= CS_COMBUSTION_LW;
1552 m_type++) {
1553 if (cs_glob_physical_model_flag[m_type] > -1)
1554 gas_combustion = true;
1555 }
1556
1557 /* First pass only: initialize izfppp */
1558
1559 static bool initialized = false;
1560 if (initialized == false) {
1561 _init_zones(n_b_faces, nozppm, izfppp);
1562 initialized = true;
1563 }
1564
1565 cs_tree_node_t *tn_bc = cs_tree_get_node(cs_glob_tree,
1566 "boundary_conditions/boundary");
1567
1568 #if _XML_DEBUG_
1569 bft_printf("==> %s\n", __func__);
1570 bft_printf("--boundary zones count: %i\n", boundaries->n_zones);
1571 #endif
1572
1573 /* At each time-step, loop on boundary face zones:
1574 set itypfb, rcodcl and icodcl thanks to the arrays
1575 of the structures defined in the first part of the function */
1576
1577 for (int izone = 0; izone < boundaries->n_zones; izone++) {
1578
1579 int zone_nbr = boundaries->bc_num[izone];
1580 const cs_zone_t *bz = cs_boundary_zone_by_id(zone_nbr);
1581
1582 #if _XML_DEBUG_
1583 bft_printf("\n---zone %i label: %s\n", zone_nbr, boundaries->label[izone]);
1584 bft_printf("---zone %i nature: %s\n", zone_nbr, boundaries->nature[izone]);
1585 bft_printf("---zone %i number of faces: %i\n", zone_nbr, bz->n_elts);
1586 #endif
1587
1588 int wall_type = 1;
1589
1590 /* Mapped inlet? */
1591
1592 if (cs_gui_strcmp(boundaries->nature[izone], "wall")) {
1593 if (boundaries->rough[izone] >= 0.0)
1594 wall_type = 6; //TODO remove and use all roughness wall function
1595 else
1596 wall_type = 5;
1597 }
1598
1599 /* for each field */
1600 for (int f_id = 0; f_id < n_fields; f_id++) {
1601 const cs_field_t *f = cs_field_by_id(f_id);
1602 const int var_key_id = cs_field_key_id("variable_id");
1603 cs_lnum_t ivar = cs_field_get_key_int(f, var_key_id) -1;
1604
1605 if (f->type & CS_FIELD_CDO)
1606 continue; /* TODO: Avoid a SIGSEV (when sharing CDO and FV, one has to
1607 find a better fix) */
1608
1609 if (f->type & CS_FIELD_VARIABLE) {
1610
1611 switch (boundaries->type_code[f->id][izone]) {
1612
1613 case DIRICHLET_CNV:
1614 {
1615 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1616 cs_lnum_t face_id = bz->elt_ids[elt_id];
1617 for (cs_lnum_t i = 0; i < f->dim; i++) {
1618 icodcl[(ivar + i) * n_b_faces + face_id] = -wall_type;
1619 rcodcl[(ivar + i) * n_b_faces + face_id]
1620 = boundaries->values[f->id][izone * f->dim + i].val1;
1621 }
1622 }
1623 }
1624 break;
1625
1626 case DIRICHLET_FORMULA:
1627 {
1628 int icodcl_m = 1;
1629 const char *f_name = f->name;
1630
1631 if (f == CS_F_(h) && boundaries->t_to_h[izone]) {
1632 icodcl_m = -1;
1633 f_name = "temperature";
1634 }
1635
1636 cs_real_t *new_vals
1637 = cs_meg_boundary_function(bz, f_name, "dirichlet_formula");
1638
1639 for (cs_lnum_t ii = 0; ii < f->dim; ii++) {
1640 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1641 cs_lnum_t face_id = bz->elt_ids[elt_id];
1642 icodcl[(ivar + ii) *n_b_faces + face_id]
1643 = wall_type * icodcl_m;
1644 rcodcl[(ivar + ii) * n_b_faces + face_id]
1645 = new_vals[ii * bz->n_elts + elt_id];
1646 }
1647 }
1648 BFT_FREE(new_vals);
1649 break;
1650 }
1651
1652 case EXCHANGE_COEFF:
1653 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1654 cs_lnum_t face_id = bz->elt_ids[elt_id];
1655 for (cs_lnum_t i = 0; i < f->dim; i++) {
1656 icodcl[(ivar + i) * n_b_faces + face_id] = 5;
1657 rcodcl[0 * n_b_faces * (*nvar) + (ivar + i) * n_b_faces + face_id]
1658 = boundaries->values[f->id][izone * f->dim + i].val1;
1659 rcodcl[1 * n_b_faces * (*nvar) + (ivar + i) * n_b_faces + face_id]
1660 = boundaries->values[f->id][izone * f->dim + i].val2;
1661 }
1662 }
1663 break;
1664
1665 case NEUMANN_FORMULA:
1666 {
1667 cs_real_t *new_vals =
1668 cs_meg_boundary_function(bz, f->name, "neumann_formula");
1669
1670 for (cs_lnum_t ii = 0; ii < f->dim; ii++) {
1671 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1672 cs_lnum_t face_id = bz->elt_ids[elt_id];
1673 icodcl[(ivar + ii) *n_b_faces + face_id] = 3;
1674 rcodcl[2 * n_b_faces * (*nvar) + (ivar + ii) * n_b_faces + face_id]
1675 = new_vals[ii * bz->n_elts + elt_id];
1676 }
1677 }
1678 BFT_FREE(new_vals);
1679 break;
1680 }
1681
1682 case EXCHANGE_COEFF_FORMULA:
1683 {
1684 cs_real_t *new_vals =
1685 cs_meg_boundary_function(bz, f->name, "exchange_coefficient_formula");
1686
1687 for (cs_lnum_t ii = 0; ii < f->dim; ii++) {
1688 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1689 cs_lnum_t face_id = bz->elt_ids[elt_id];
1690 icodcl[(ivar + ii) *n_b_faces + face_id] = 5;
1691
1692 rcodcl[0 * n_b_faces * (*nvar) + (ivar + ii) * n_b_faces + face_id]
1693 = new_vals[ii * bz->n_elts + elt_id];
1694
1695 rcodcl[1 * n_b_faces * (*nvar) + (ivar + ii) * n_b_faces + face_id]
1696 = new_vals[f->dim * bz->n_elts + elt_id];
1697 }
1698 }
1699 BFT_FREE(new_vals);
1700 break;
1701 }
1702 }
1703 } /* switch */
1704 } /* Loop on fields */
1705
1706 if (cs_glob_physical_model_flag[CS_JOULE_EFFECT] > -1) {
1707 if (cs_glob_elec_option->ielcor == 1) {
1708 const cs_field_t *f = CS_F_(potr);
1709 const int var_key_id = cs_field_key_id("variable_id");
1710 cs_lnum_t ivar = cs_field_get_key_int(f, var_key_id) -1;
1711
1712 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1713 cs_lnum_t face_id = bz->elt_ids[elt_id];
1714 rcodcl[ivar * n_b_faces + face_id] *= cs_glob_elec_option->coejou;
1715 }
1716
1717 int ieljou = cs_glob_physical_model_flag[CS_JOULE_EFFECT];
1718 if (ieljou == 2 || ieljou == 4) {
1719 const cs_field_t *fi = CS_F_(poti);
1720 ivar = cs_field_get_key_int(fi, var_key_id) -1;
1721 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1722 cs_lnum_t face_id = bz->elt_ids[elt_id];
1723 rcodcl[ivar * n_b_faces + face_id] *= cs_glob_elec_option->coejou;
1724 }
1725 }
1726 }
1727 }
1728
1729 if (cs_glob_physical_model_flag[CS_ELECTRIC_ARCS] > -1) {
1730 const cs_field_t *f = CS_F_(potr);
1731 const int var_key_id = cs_field_key_id("variable_id");
1732 cs_lnum_t ivar = cs_field_get_key_int(f, var_key_id) -1;
1733
1734 if ( boundaries->type_code[f->id][izone] == DIRICHLET_IMPLICIT
1735 && cs_glob_elec_option->ielcor == 1) {
1736 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1737 cs_lnum_t face_id = bz->elt_ids[elt_id];
1738 icodcl[ivar * n_b_faces + face_id] = 5;
1739 rcodcl[ivar * n_b_faces + face_id] = cs_glob_elec_option->pot_diff;
1740 }
1741 }
1742
1743 const cs_field_t *fp = cs_field_by_name_try("vec_potential");
1744 ivar = cs_field_get_key_int(fp, var_key_id) -1;
1745
1746 if (boundaries->type_code[fp->id][izone] == NEUMANN_IMPLICIT)
1747 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1748 cs_lnum_t face_id = bz->elt_ids[elt_id];
1749 cs_lnum_t iel = b_face_cells[face_id];
1750 icodcl[ivar *n_b_faces + face_id] = 5;
1751 icodcl[(ivar+1) *n_b_faces + face_id] = 5;
1752 icodcl[(ivar+2) *n_b_faces + face_id] = 5;
1753 rcodcl[ivar * n_b_faces + face_id] = fp->val_pre[3*iel];
1754 rcodcl[(ivar+1) * n_b_faces + face_id] = fp->val_pre[3*iel+1];
1755 rcodcl[(ivar+2) * n_b_faces + face_id] = fp->val_pre[3*iel+2];
1756 }
1757 }
1758
1759 /* Boundary conditions by boundary type
1760 ------------------------------------ */
1761
1762 if (cs_gui_strcmp(boundaries->nature[izone], "inlet")) {
1763
1764 tn_bc = _get_zone_bc_node(tn_bc, izone);
1765
1766 cs_tree_node_t *tn_vp = cs_tree_node_get_child(tn_bc, "velocity_pressure");
1767 const char *choice_v = cs_gui_node_get_tag(tn_vp, "choice");
1768 const char *choice_d = cs_gui_node_get_tag(tn_vp, "direction");
1769
1770 /* Update the zone's arrays (iqimp, dh, xintur, icalke, qimp,...)
1771 because they are re-initialized at each time step
1772 in PRECLI and PPPRCL routines */
1773
1774 /* data by zone */
1775 iqimp[zone_nbr-1] = boundaries->iqimp[izone];
1776 dh[zone_nbr-1] = boundaries->dh[izone];
1777 xintur[zone_nbr-1] = boundaries->xintur[izone];
1778 icalke[zone_nbr-1] = boundaries->icalke[izone];
1779
1780 if (solid_fuels) {
1781 const cs_combustion_model_t *cm = cs_glob_combustion_model;
1782
1783 ientat[zone_nbr-1] = boundaries->ientat[izone];
1784 inmoxy[zone_nbr-1] = boundaries->inmoxy[izone];
1785 ientcp[zone_nbr-1] = boundaries->ientcp[izone];
1786 qimpat[zone_nbr-1] = boundaries->qimp[izone];
1787 timpat[zone_nbr-1] = boundaries->timpat[izone];
1788
1789 for (int icharb = 0; icharb < cm->coal.n_coals; icharb++) {
1790 int ich = icharb *(*nozppm)+zone_nbr-1;
1791 qimpcp[ich] = boundaries->qimpcp[izone][icharb];
1792 timpcp[ich] = boundaries->timpcp[izone][icharb];
1793
1794 for (int iclass = 0;
1795 iclass < cm->coal.n_classes_per_coal[icharb];
1796 iclass++) {
1797 int icl = iclass * (*nozppm) * ncharm + ich;
1798 distch[icl] = boundaries->distch[izone][icharb][iclass];
1799 }
1800 }
1801 }
1802 else if (gas_combustion) {
1803 ientfu[zone_nbr-1] = boundaries->ientfu[izone];
1804 ientox[zone_nbr-1] = boundaries->ientox[izone];
1805 ientgb[zone_nbr-1] = boundaries->ientgb[izone];
1806 ientgf[zone_nbr-1] = boundaries->ientgf[izone];
1807 tkent[zone_nbr-1] = boundaries->tkent[izone];
1808 fment[zone_nbr-1] = boundaries->fment[izone];
1809
1810 if ( cs_gui_strcmp(choice_v, "flow1_formula")
1811 || cs_gui_strcmp(choice_v, "flow2_formula")) {
1812
1813 if (cs_gui_strcmp(choice_v, "flow1_formula")) {
1814 qimp[zone_nbr-1] = *cs_meg_boundary_function(bz,
1815 "velocity",
1816 "flow1_formula");
1817 } else if (cs_gui_strcmp(choice_v, "flow2_formula")) {
1818 qimp[zone_nbr-1] = *cs_meg_boundary_function(bz,
1819 "velocity",
1820 "flow2_formula");
1821 }
1822 }
1823 else {
1824 qimp[zone_nbr-1] = boundaries->qimp[izone];
1825 }
1826 }
1827 else if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
1828 const int var_key_id = cs_field_key_id("variable_id");
1829
1830 if ( boundaries->itype[izone] == CS_ESICF
1831 ||boundaries->itype[izone] == CS_EPHCF) {
1832 const cs_field_t *fp = cs_field_by_name_try("pressure");
1833 int ivarp = cs_field_get_key_int(fp, var_key_id) -1;
1834
1835 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1836 cs_lnum_t face_id = bz->elt_ids[elt_id];
1837 rcodcl[ivarp * n_b_faces + face_id] = boundaries->prein[izone];
1838 }
1839 qimp[zone_nbr-1] = boundaries->qimp[izone];
1840 }
1841
1842 if (boundaries->itype[izone] == CS_ESICF) {
1843 cs_field_t *b_rho = cs_field_by_name_try("boundary_density");
1844 const cs_field_t *ft = cs_field_by_name_try("temperature");
1845 int ivart = cs_field_get_key_int(ft, var_key_id) -1;
1846
1847 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1848 cs_lnum_t face_id = bz->elt_ids[elt_id];
1849 rcodcl[ivart * n_b_faces + face_id] = boundaries->tempin[izone];
1850 b_rho->val[face_id] = boundaries->rhoin[izone];
1851 }
1852 }
1853 }
1854 else {
1855 if (boundaries->velocity_e[izone]) {
1856 if (cs_gui_strcmp(choice_v, "flow1_formula")) {
1857 qimp[zone_nbr-1] = *cs_meg_boundary_function(bz,
1858 "velocity",
1859 "flow1_formula");
1860 }
1861 else if (cs_gui_strcmp(choice_v, "flow2_formula")) {
1862 qimp[zone_nbr-1] = *cs_meg_boundary_function(bz,
1863 "velocity",
1864 "flow2_formula");
1865 }
1866 }
1867 else {
1868 qimp[zone_nbr-1] = boundaries->qimp[izone];
1869 }
1870 }
1871
1872 /* data by boundary faces */
1873
1874 int inlet_type = CS_INLET;
1875
1876 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1)
1877 inlet_type = boundaries->itype[izone];
1878 else {
1879 int convective_inlet = 0;
1880 _boundary_status("inlet", boundaries->label[izone],
1881 "convective_inlet", &convective_inlet);
1882 if (convective_inlet)
1883 inlet_type = CS_CONVECTIVE_INLET;
1884 }
1885
1886 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1887 cs_lnum_t face_id = bz->elt_ids[elt_id];
1888
1889 /* zone number and nature of boundary */
1890 izfppp[face_id] = zone_nbr;
1891 itypfb[face_id] = inlet_type;
1892 }
1893
1894 if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] > -1) {
1895 iprofm[zone_nbr-1] = boundaries->meteo[izone].read_data;
1896 if (iprofm[zone_nbr-1] == 1) {
1897 choice_v = NULL;
1898 choice_d = NULL;
1899 }
1900 if (boundaries->meteo[izone].automatic) {
1901 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1902 cs_lnum_t face_id = bz->elt_ids[elt_id];
1903 iautom[face_id] = 1;
1904 }
1905 }
1906 }
1907
1908 /* Dirichlet for velocity */
1909
1910 const cs_field_t *fv = cs_field_by_name_try("velocity");
1911 const int var_key_id = cs_field_key_id("variable_id");
1912 int ivarv = cs_field_get_key_int(fv, var_key_id) -1;
1913 if (cs_gui_strcmp(choice_d, "coordinates")) {
1914 if (cs_gui_strcmp(choice_v, "norm")) {
1915 norm = boundaries->norm[izone]
1916 / cs_math_3_norm(boundaries->dir[izone]);
1917
1918 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1919 cs_lnum_t face_id = bz->elt_ids[elt_id];
1920 for (cs_lnum_t ic = 0; ic < 3; ic++) {
1921 rcodcl[(ivarv + ic) * n_b_faces + face_id]
1922 = boundaries->dir[izone][ic] * norm;
1923 }
1924 }
1925 }
1926 else if ( cs_gui_strcmp(choice_v, "flow1")
1927 || cs_gui_strcmp(choice_v, "flow2")
1928 || cs_gui_strcmp(choice_v, "flow1_formula")
1929 || cs_gui_strcmp(choice_v, "flow2_formula")) {
1930 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1931 cs_lnum_t face_id = bz->elt_ids[elt_id];
1932 for (cs_lnum_t ic = 0; ic < 3; ic++) {
1933 rcodcl[(ivarv + ic) * n_b_faces + face_id]
1934 = boundaries->dir[izone][ic];
1935 }
1936 }
1937 }
1938 else if (cs_gui_strcmp(choice_v, "norm_formula")) {
1939 cs_real_t *new_vals = cs_meg_boundary_function(bz,
1940 "velocity",
1941 "norm_formula");
1942
1943 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1944 cs_lnum_t face_id = bz->elt_ids[elt_id];
1945
1946 cs_real_t x_norm = cs_math_3_norm(boundaries->dir[izone]);
1947 if (x_norm <= 0.)
1948 bft_error(__FILE__, __LINE__, 0,
1949 _("Error in the boundary conditions: "
1950 " the normal direction is of norm 0."));
1951
1952 for (cs_lnum_t ic = 0; ic < 3; ic++) {
1953 rcodcl[(ivarv + ic) * n_b_faces + face_id]
1954 = boundaries->dir[izone][ic] *
1955 new_vals[elt_id] / x_norm;
1956 }
1957 }
1958 BFT_FREE(new_vals);
1959 }
1960 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
1961 if (boundaries->itype[izone] == CS_EPHCF) {
1962 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1963 cs_lnum_t face_id = bz->elt_ids[elt_id];
1964 for (cs_lnum_t ic = 0; ic < 3; ic++) {
1965 rcodcl[(ivarv + ic) * n_b_faces + face_id]
1966 = boundaries->dir[izone][ic];
1967 }
1968 }
1969 }
1970 }
1971 }
1972 else if ( cs_gui_strcmp(choice_d, "normal")
1973 || cs_gui_strcmp(choice_d, "translation")) {
1974 if (cs_gui_strcmp(choice_v, "norm")) {
1975 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1976 cs_lnum_t face_id = bz->elt_ids[elt_id];
1977
1978 norm = boundaries->norm[izone] / b_face_surf[face_id];
1979
1980 for (cs_lnum_t i = 0; i < 3; i++)
1981 rcodcl[(ivarv + i) * n_b_faces + face_id]
1982 = -b_face_normal[face_id][i] * norm;
1983 }
1984 }
1985 else if ( cs_gui_strcmp(choice_v, "flow1")
1986 || cs_gui_strcmp(choice_v, "flow2")
1987 || cs_gui_strcmp(choice_v, "flow1_formula")
1988 || cs_gui_strcmp(choice_v, "flow2_formula")) {
1989 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
1990 cs_lnum_t face_id = bz->elt_ids[elt_id];
1991
1992 for (cs_lnum_t i = 0; i < 3; i++)
1993 rcodcl[(ivarv + i) * n_b_faces + face_id]
1994 = -b_face_normal[face_id][i] / b_face_surf[face_id];
1995 }
1996 }
1997 else if (cs_gui_strcmp(choice_v, "norm_formula")) {
1998 cs_real_t *new_vals = cs_meg_boundary_function(bz,
1999 "velocity",
2000 "norm_formula");
2001
2002 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2003 cs_lnum_t face_id = bz->elt_ids[elt_id];
2004
2005 for (cs_lnum_t i = 0; i < 3; i++)
2006 rcodcl[(ivarv + i) * n_b_faces + face_id]
2007 = -b_face_normal[face_id][i] *
2008 new_vals[elt_id] / b_face_surf[face_id];
2009 }
2010 BFT_FREE(new_vals);
2011 }
2012
2013 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
2014 if (boundaries->itype[izone] == CS_EPHCF) {
2015 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2016 cs_lnum_t face_id = bz->elt_ids[elt_id];
2017
2018 for (cs_lnum_t i = 0; i < 3; i++)
2019 rcodcl[(ivarv + i) * n_b_faces + face_id]
2020 = -b_face_normal[face_id][i];
2021 }
2022 }
2023 }
2024 }
2025 else if (cs_gui_strcmp(choice_d, "formula")) {
2026 cs_real_t *xvals = cs_meg_boundary_function(bz,
2027 "direction",
2028 "formula");
2029
2030 if (cs_gui_strcmp(choice_v, "norm")) {
2031 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2032 cs_lnum_t face_id = bz->elt_ids[elt_id];
2033
2034 cs_real_t x[3];
2035 for (int ic = 0; ic < 3; ic++)
2036 x[ic] = xvals[ic*bz->n_elts + elt_id];
2037
2038 cs_real_t x_norm = cs_math_3_norm(x);
2039
2040 if (x_norm <= 0.)
2041 bft_error(__FILE__, __LINE__, 0,
2042 _("Error in the boundary conditions: "
2043 "the normal direction is of norm 0.\n "));
2044
2045 norm = boundaries->norm[izone] / x_norm;
2046
2047 for (int i = 0; i < 3; i++)
2048 rcodcl[(ivarv + i) * n_b_faces + face_id] = x[i] * norm;
2049 }
2050 }
2051 else if ( cs_gui_strcmp(choice_v, "flow1")
2052 || cs_gui_strcmp(choice_v, "flow2")
2053 || cs_gui_strcmp(choice_v, "flow1_formula")
2054 || cs_gui_strcmp(choice_v, "flow2_formula")) {
2055 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2056 cs_lnum_t face_id = bz->elt_ids[elt_id];
2057
2058 cs_real_t x[3];
2059 for (int ic = 0; ic < 3; ic++)
2060 x[ic] = xvals[ic*bz->n_elts + elt_id];
2061
2062 cs_real_t x_norm = cs_math_3_norm(x);
2063 if (x_norm <= 0.)
2064 bft_error(__FILE__, __LINE__, 0,
2065 _("Error in the boundary conditions: "
2066 "the normal direction is of norm 0.\n "));
2067
2068 for (int i = 0; i < 3; i++)
2069 rcodcl[(ivarv + i) * n_b_faces + face_id] = x[i];
2070 }
2071 }
2072 else if (cs_gui_strcmp(choice_v, "norm_formula")) {
2073 cs_real_t *norm_vals = cs_meg_boundary_function(bz,
2074 "velocity",
2075 "norm_formula");
2076
2077 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2078 cs_lnum_t face_id = bz->elt_ids[elt_id];
2079
2080 cs_real_t x[3];
2081 for (int ic = 0; ic < 3; ic++)
2082 x[ic] = xvals[ic*bz->n_elts + elt_id];
2083
2084 cs_real_t x_norm = cs_math_3_norm(x);
2085 if (x_norm <= 0.)
2086 bft_error(__FILE__, __LINE__, 0,
2087 _("Error in the boundary conditions: "
2088 "the normal direction is of norm 0.\n "));
2089
2090 norm = norm_vals[elt_id] / x_norm;
2091
2092 for (int i = 0; i < 3; i++)
2093 rcodcl[(ivarv + i) * n_b_faces + face_id] = x[i] * norm;
2094 }
2095 }
2096
2097 BFT_FREE(xvals);
2098
2099 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
2100 if (boundaries->itype[izone] == CS_EPHCF) {
2101 xvals = cs_meg_boundary_function(bz,
2102 "direction",
2103 "formula");
2104 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2105 cs_lnum_t face_id = bz->elt_ids[elt_id];
2106
2107 cs_real_t x[3];
2108 for (int ic = 0; ic < 3; ic++)
2109 x[ic] = xvals[ic*bz->n_elts + elt_id];
2110
2111 for (cs_lnum_t i = 0; i < 3; i++)
2112 rcodcl[(ivarv + i) * n_b_faces + face_id] = x[i];
2113 }
2114 BFT_FREE(xvals);
2115 }
2116 }
2117 }
2118
2119 /* turbulent inlet, with formula */
2120 if (icalke[zone_nbr-1] == 0) {
2121
2122 tn_bc = _get_zone_bc_node(tn_bc, izone);
2123
2124 cs_tree_node_t *tn_t = cs_tree_node_get_child(tn_bc, "turbulence");
2125 const char *formula = cs_tree_node_get_child_value_str(tn_t, "formula");
2126
2127 if (formula != NULL) {
2128
2129 const char *model = cs_gui_get_thermophysical_model("turbulence");
2130 if (model == NULL)
2131 return;
2132
2133 if ( cs_gui_strcmp(model, "k-epsilon")
2134 || cs_gui_strcmp(model, "k-epsilon-PL")) {
2135
2136 cs_real_t *new_vals = cs_meg_boundary_function(bz,
2137 "turbulence_ke",
2138 "formula");
2139
2140 cs_field_t *c_k = cs_field_by_name("k");
2141 cs_field_t *c_eps = cs_field_by_name("epsilon");
2142 int ivark = cs_field_get_key_int(c_k, var_key_id) -1;
2143 int ivare = cs_field_get_key_int(c_eps, var_key_id) -1;
2144
2145 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2146 cs_lnum_t face_id = bz->elt_ids[elt_id];
2147 rcodcl[ivark * n_b_faces + face_id]
2148 = new_vals[0 * bz->n_elts + elt_id];
2149 rcodcl[ivare * n_b_faces + face_id]
2150 = new_vals[1 * bz->n_elts + elt_id];
2151 }
2152 BFT_FREE(new_vals);
2153 }
2154 else if ( cs_gui_strcmp(model, "Rij-epsilon")
2155 ||cs_gui_strcmp(model, "Rij-SSG")) {
2156
2157 cs_real_t *new_vals = cs_meg_boundary_function(bz,
2158 "turbulence_rije",
2159 "formula");
2160 cs_field_t *cfld_rij = cs_field_by_name("rij");
2161 cs_field_t *c_eps = cs_field_by_name("epsilon");
2162 int ivarrij = cs_field_get_key_int(cfld_rij, var_key_id) - 1;
2163 int ivare = cs_field_get_key_int(c_eps, var_key_id) -1;
2164
2165 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2166 cs_lnum_t face_id = bz->elt_ids[elt_id];
2167
2168 /* Values are stored for rij components then epsilon */
2169 for (int ii = 0; ii < 6; ii++)
2170 rcodcl[(ivarrij + ii) * n_b_faces + face_id]
2171 = new_vals[bz->n_elts * ii + elt_id];
2172
2173 rcodcl[ivare * n_b_faces + face_id]
2174 = new_vals[bz->n_elts * 6 + elt_id];
2175 }
2176 BFT_FREE(new_vals);
2177 }
2178 else if (cs_gui_strcmp(model, "Rij-EBRSM")) {
2179
2180 cs_real_t *new_vals = cs_meg_boundary_function(bz,
2181 "turbulence_rij_ebrsm",
2182 "formula");
2183
2184 cs_field_t *cfld_rij = cs_field_by_name("rij");
2185 cs_field_t *c_eps = cs_field_by_name("epsilon");
2186 cs_field_t *c_a = cs_field_by_name("alpha");
2187
2188 int ivarrij = cs_field_get_key_int(cfld_rij, var_key_id) - 1;
2189 int ivare = cs_field_get_key_int(c_eps, var_key_id) -1;
2190 int ivara = cs_field_get_key_int(c_a, var_key_id) -1;
2191
2192 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2193 cs_lnum_t face_id = bz->elt_ids[elt_id];
2194
2195 /* Values are stored for rij components then epsilon and alpha*/
2196 for (int ii = 0; ii < 6; ii++)
2197 rcodcl[(ivarrij + ii) * n_b_faces + face_id]
2198 = new_vals[bz->n_elts * ii + elt_id];
2199
2200 rcodcl[ivare * n_b_faces + face_id]
2201 = new_vals[bz->n_elts * 6 + elt_id];
2202
2203 rcodcl[ivara * n_b_faces + face_id]
2204 = new_vals[bz->n_elts * 7 + elt_id];
2205 }
2206 BFT_FREE(new_vals);
2207 }
2208 else if (cs_gui_strcmp(model, "v2f-BL-v2/k")) {
2209 cs_real_t *new_vals = cs_meg_boundary_function(bz,
2210 "turbulence_v2f",
2211 "formula");
2212
2213 cs_field_t *c_k = cs_field_by_name("k");
2214 cs_field_t *c_eps = cs_field_by_name("epsilon");
2215 cs_field_t *c_phi = cs_field_by_name("phi");
2216 cs_field_t *c_a = cs_field_by_name("alpha");
2217 int ivark = cs_field_get_key_int(c_k, var_key_id) -1;
2218 int ivare = cs_field_get_key_int(c_eps, var_key_id) -1;
2219 int ivarp = cs_field_get_key_int(c_phi, var_key_id) -1;
2220 int ivara = cs_field_get_key_int(c_a, var_key_id) -1;
2221
2222 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2223 cs_lnum_t face_id = bz->elt_ids[elt_id];
2224 rcodcl[ivark * n_b_faces + face_id]
2225 = new_vals[0 * bz->n_elts + elt_id];
2226 rcodcl[ivare * n_b_faces + face_id]
2227 = new_vals[1 * bz->n_elts + elt_id];
2228 rcodcl[ivarp * n_b_faces + face_id]
2229 = new_vals[2 * bz->n_elts + elt_id];
2230 rcodcl[ivara * n_b_faces + face_id]
2231 = new_vals[3 * bz->n_elts + elt_id];
2232 }
2233 BFT_FREE(new_vals);
2234 }
2235 else if (cs_gui_strcmp(model, "k-omega-SST")) {
2236 cs_real_t *new_vals = cs_meg_boundary_function(bz,
2237 "turbulence_kw",
2238 "formula");
2239
2240 cs_field_t *c_k = cs_field_by_name("k");
2241 cs_field_t *c_o = cs_field_by_name("omega");
2242 int ivark = cs_field_get_key_int(c_k, var_key_id) -1;
2243 int ivaro = cs_field_get_key_int(c_o, var_key_id) -1;
2244
2245 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2246 cs_lnum_t face_id = bz->elt_ids[elt_id];
2247 rcodcl[ivark * n_b_faces + face_id]
2248 = new_vals[0 * bz->n_elts + elt_id];
2249 rcodcl[ivaro * n_b_faces + face_id]
2250 = new_vals[1 * bz->n_elts + elt_id];
2251 }
2252 BFT_FREE(new_vals);
2253 }
2254 else if (cs_gui_strcmp(model, "Spalart-Allmaras")) {
2255 cs_real_t *new_vals = cs_meg_boundary_function(bz,
2256 "turbulence_spalart",
2257 "formula");
2258
2259 cs_field_t *c_nu = cs_field_by_name("nu_tilda");
2260 int ivarnu = cs_field_get_key_int(c_nu, var_key_id) -1;
2261
2262 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2263 cs_lnum_t face_id = bz->elt_ids[elt_id];
2264 rcodcl[ivarnu * n_b_faces + face_id] = new_vals[elt_id];
2265 }
2266 BFT_FREE(new_vals);
2267 }
2268 else
2269 bft_error(__FILE__, __LINE__, 0,
2270 _("Invalid turbulence model: %s.\n"), model);
2271 }
2272 }
2273
2274 #if _XML_DEBUG_
2275 if (cs_gui_strcmp(choice_v, "norm"))
2276 bft_printf("-----velocity: %s => %12.5e \n",
2277 choice_v, boundaries->norm[izone]);
2278 if (cs_gui_strcmp(choice_v, "flow1") || cs_gui_strcmp(choice_v, "flow2"))
2279 bft_printf("-----velocity: %s => %12.5e \n",
2280 choice_v, boundaries->qimp[izone]);
2281 if ( cs_gui_strcmp(choice_v, "norm_formula")
2282 || cs_gui_strcmp(choice_v, "flow1_formula")
2283 || cs_gui_strcmp(choice_v, "flow2_formula"))
2284 bft_printf("-----velocity: %s => %d \n",
2285 choice_v, (boundaries->velocity_e[izone] ? 1: 0));
2286 if ( cs_gui_strcmp(choice_d, "coordinates")
2287 || cs_gui_strcmp(choice_d, "translation"))
2288 bft_printf("-----direction: %s => %12.5e %12.5e %12.5e\n",
2289 choice_v, boundaries->dir[izone][0],
2290 boundaries->dir[izone][1], boundaries->dir[izone][2]);
2291
2292 if (solid_fuels) {
2293 bft_printf("-----iqimp=%i, qimpat=%12.5e \n",
2294 iqimp[zone_nbr-1], qimpat[zone_nbr-1]);
2295 bft_printf("-----ientat=%i, ientcp=%i, timpat=%12.5e \n",
2296 ientat[zone_nbr-1], ientcp[zone_nbr-1], timpat[zone_nbr-1]);
2297
2298 for (int icharb = 0; icharb < boundaries->n_coals; icharb++) {
2299 bft_printf("-----coal=%i, qimpcp=%12.5e, timpcp=%12.5e \n",
2300 icharb+1, qimpcp[icharb *(*nozppm)+zone_nbr-1],
2301 timpcp[icharb *(*nozppm)+zone_nbr-1]);
2302
2303 for (int iclass = 0; iclass < nclpch[icharb]; iclass++)
2304 bft_printf("-----coal=%i, class=%i, distch=%f \n",
2305 icharb+1, iclass+1,
2306 distch[ iclass * (*nozppm) * ncharm
2307 + icharb * (*nozppm) +zone_nbr-1]);
2308 }
2309 }
2310 else if (gas_combustion) {
2311 bft_printf("-----iqimp=%i \n",
2312 iqimp[zone_nbr-1]);
2313 bft_printf("-----ientox=%i, ientfu=%i, ientgf=%i, ientgb=%i \n",
2314 ientox[zone_nbr-1], ientfu[zone_nbr-1],
2315 ientgf[zone_nbr-1], ientgb[zone_nbr-1]);
2316 }
2317 else if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
2318 if (boundaries->itype[izone] == CS_ESICF) {
2319 bft_printf("-----imposed_inlet\n");
2320 bft_printf("-----premin=%g \n",boundaries->prein[zone_nbr-1]);
2321 bft_printf("-----rhoin=%g \n",boundaries->rhoin[zone_nbr-1]);
2322 bft_printf("-----tempin=%g \n",boundaries->tempin[zone_nbr-1]);
2323 }
2324 if (boundaries->itype[izone] == CS_EPHCF) {
2325 bft_printf("-----subsonic_inlet_PH\n");
2326 bft_printf("-----prein=%g \n",boundaries->prein[zone_nbr-1]);
2327 }
2328 }
2329 else {
2330 bft_printf("-----iqimp=%i, qimp=%12.5e \n",
2331 iqimp[zone_nbr-1], qimp[zone_nbr-1]);
2332 }
2333 bft_printf("-----icalke=%i, dh=%12.5e, xintur=%12.5e \n",
2334 icalke[zone_nbr-1], dh[zone_nbr-1], xintur[zone_nbr-1]);
2335
2336 if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] > -1) {
2337 bft_printf("-----iprofm=%i, automatic=%i \n",
2338 iprofm[zone_nbr-1], boundaries->meteo[izone].automatic);
2339 #endif
2340
2341 }
2342
2343 else if (cs_gui_strcmp(boundaries->nature[izone], "wall")) {
2344 int iwall = CS_SMOOTHWALL;
2345
2346 if (boundaries->rough[izone] >= 0.0) {
2347
2348 iwall = CS_ROUGHWALL;
2349 cs_field_t *f_roughness = cs_field_by_name_try("boundary_roughness");
2350 cs_field_t *f_roughness_t
2351 = cs_field_by_name_try("boundary_thermal_roughness");
2352
2353 /* Roughness value (z0) */
2354 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2355 cs_lnum_t face_id = bz->elt_ids[elt_id];
2356 assert(f_roughness != NULL);
2357 f_roughness->val[face_id] = boundaries->rough[izone];
2358
2359 /* Thermal Roughness value.
2360 * In this case thermal roughness is equal to the roughness. */
2361 if (f_roughness_t != NULL)
2362 f_roughness_t->val[face_id] = boundaries->rough[izone];
2363 }
2364 }
2365
2366 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2367 cs_lnum_t face_id = bz->elt_ids[elt_id];
2368 izfppp[face_id] = zone_nbr;
2369 itypfb[face_id] = iwall;
2370 }
2371 }
2372
2373 else if (cs_gui_strcmp(boundaries->nature[izone], "outlet")) {
2374 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2375 cs_lnum_t face_id = bz->elt_ids[elt_id];
2376 izfppp[face_id] = zone_nbr;
2377 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1)
2378 itypfb[face_id] = boundaries->itype[izone];
2379 else
2380 itypfb[face_id] = CS_OUTLET;
2381 }
2382
2383 if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] > -1) {
2384 iprofm[zone_nbr-1] = boundaries->meteo[izone].read_data;
2385 if (boundaries->meteo[izone].automatic) {
2386 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2387 cs_lnum_t face_id = bz->elt_ids[elt_id];
2388 iautom[face_id] = 1;
2389 }
2390 }
2391 }
2392 }
2393
2394 else if (cs_gui_strcmp(boundaries->nature[izone], "imposed_p_outlet")) {
2395 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2396 cs_lnum_t face_id = bz->elt_ids[elt_id];
2397 izfppp[face_id] = zone_nbr;
2398 itypfb[face_id] = CS_OUTLET;
2399 }
2400 }
2401
2402 else if (cs_gui_strcmp(boundaries->nature[izone], "symmetry")) {
2403 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2404 cs_lnum_t face_id = bz->elt_ids[elt_id];
2405 izfppp[face_id] = zone_nbr;
2406 itypfb[face_id] = CS_SYMMETRY;
2407 }
2408 }
2409
2410 else if (cs_gui_strcmp(boundaries->nature[izone], "free_inlet_outlet")) {
2411 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2412 cs_lnum_t face_id = bz->elt_ids[elt_id];
2413 izfppp[face_id] = zone_nbr;
2414 itypfb[face_id] = CS_FREE_INLET;
2415 }
2416
2417 if (boundaries->head_loss_e[izone]) {
2418 cs_real_t *new_vals =
2419 cs_meg_boundary_function(bz, "head_loss", "formula");
2420
2421 const cs_field_t *fp = cs_field_by_name_try("pressure");
2422 const int var_key_id = cs_field_key_id("variable_id");
2423 int ivarp = cs_field_get_key_int(fp, var_key_id) -1;
2424
2425 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2426 cs_lnum_t face_id = bz->elt_ids[elt_id];
2427 rcodcl[1 * n_b_faces * (*nvar) + ivarp * n_b_faces + face_id]
2428 = new_vals[elt_id];
2429 }
2430 BFT_FREE(new_vals);
2431 }
2432 }
2433
2434 else if (cs_gui_strcmp(boundaries->nature[izone], "free_surface")) {
2435 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2436 cs_lnum_t face_id = bz->elt_ids[elt_id];
2437 izfppp[face_id] = zone_nbr;
2438 itypfb[face_id] = CS_FREE_SURFACE;
2439 }
2440 }
2441
2442 else if (cs_gui_strcmp(boundaries->nature[izone], "groundwater")) {
2443
2444 const int var_key_id = cs_field_key_id("variable_id");
2445
2446 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2447 cs_lnum_t face_id = bz->elt_ids[elt_id];
2448 izfppp[face_id] = zone_nbr;
2449 itypfb[face_id] = CS_INDEF;
2450 }
2451
2452 int ivar1 = -1;
2453 const cs_field_t *fp1 = cs_field_by_name_try("hydraulic_head");
2454 if (fp1 != NULL)
2455 ivar1 = cs_field_get_key_int(fp1, var_key_id) -1;
2456
2457 /* set velocity to 0 */
2458 const cs_field_t *fp2 = cs_field_by_name_try("velocity");
2459 if (fp2 != NULL) {
2460 int ivar2 = cs_field_get_key_int(fp2, var_key_id) -1;
2461
2462 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2463 cs_lnum_t face_id = bz->elt_ids[elt_id];
2464 for (cs_lnum_t i = 0; i < 3; i++) {
2465 icodcl[(ivar2 + i) * n_b_faces + face_id] = 3;
2466 rcodcl[(ivar2 + i) * n_b_faces + face_id] = 0.;
2467 }
2468 }
2469 }
2470
2471 if (ivar1 > -1) { /* groundwater model is active */
2472
2473 tn_bc = _get_zone_bc_node(tn_bc, izone);
2474 cs_tree_node_t *tn_hh = cs_tree_node_get_child(tn_bc, "hydraulicHead");
2475 const char *choice_d = cs_gui_node_get_tag(tn_hh, "choice");
2476
2477 if (cs_gui_strcmp(choice_d, "dirichlet_formula")) {
2478 cs_real_t *new_vals = cs_meg_boundary_function(bz,
2479 "hydraulic_head",
2480 "dirichlet_formula");
2481
2482 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2483 cs_lnum_t face_id = bz->elt_ids[elt_id];
2484 icodcl[ivar1 * n_b_faces + face_id] = 1;
2485 rcodcl[ivar1 * n_b_faces + face_id] = new_vals[elt_id];
2486 }
2487 BFT_FREE(new_vals);
2488 }
2489 }
2490
2491 }
2492
2493 else if (cs_gui_strcmp(boundaries->nature[izone], "undefined")) {
2494 for (cs_lnum_t elt_id = 0; elt_id < bz->n_elts; elt_id++) {
2495 cs_lnum_t face_id = bz->elt_ids[elt_id];
2496 izfppp[face_id] = zone_nbr;
2497 itypfb[face_id] = CS_INDEF;
2498 }
2499
2500 }
2501
2502 else {
2503 bft_error(__FILE__, __LINE__, 0,
2504 _("boundary nature %s is unknown \n"),
2505 boundaries->nature[izone]);
2506 }
2507
2508 #if _XML_DEBUG_
2509 if (bz->n_elts > 0) {
2510 cs_lnum_t face_id = bz->elt_ids[0];
2511
2512 for (int f_id = 0; f_id < n_fields; f_id++) {
2513 const cs_field_t *f = cs_field_by_id(f_id);
2514 const int var_key_id = cs_field_key_id("variable_id");
2515 cs_lnum_t ivar = cs_field_get_key_int(f, var_key_id) -1;
2516 if (f->type & CS_FIELD_VARIABLE) {
2517 bft_printf("------%s: icodcl=%i, "
2518 "rcodcl(1)=%12.5e, rcodcl(2)=%12.5e, rcodcl(3)=%12.5e\n",
2519 f->name,
2520 icodcl[ivar *n_b_faces +face_id ],
2521 rcodcl[0 * n_b_faces * (*nvar) +ivar * n_b_faces +face_id],
2522 rcodcl[1 * n_b_faces * (*nvar) +ivar * n_b_faces +face_id],
2523 rcodcl[2 * n_b_faces * (*nvar) +ivar * n_b_faces +face_id]);
2524 }
2525 }
2526 }
2527 #endif
2528
2529 } /* for (izone=0; izone < boundaries->n_zones; izone++) */
2530
2531 /* Define boundary conditions based on cs_equation_param_t structures */
2532
2533 cs_boundary_conditions_compute(*nvar,
2534 itypfb,
2535 icodcl,
2536 rcodcl);
2537 }
2538
2539 /*----------------------------------------------------------------------------
2540 * Boundary conditions input verification
2541 *
2542 * Fortran Interface:
2543 *
2544 * SUBROUTINE UICLVE
2545 * *****************
2546 *
2547 * integer nozppm <-- max number of boundary conditions zone
2548 *----------------------------------------------------------------------------*/
2549
2550 void CS_PROCF (uiclve, UICLVE)(const int *nozppm)
2551 {
2552 int inature = -1;
2553
2554 for (int izone = 0; izone < boundaries->n_zones; izone++) {
2555 if (cs_gui_strcmp(boundaries->nature[izone], "inlet")) {
2556 inature = CS_INLET;
2557 }
2558 else if (cs_gui_strcmp(boundaries->nature[izone], "wall")) {
2559 inature = CS_ROUGHWALL;
2560 if (boundaries->rough[izone] < 0.0)
2561 inature = CS_SMOOTHWALL;
2562 }
2563 else if (cs_gui_strcmp(boundaries->nature[izone], "outlet")
2564 || cs_gui_strcmp(boundaries->nature[izone], "imposed_p_outlet")) {
2565 inature = CS_OUTLET;
2566 }
2567 else if (cs_gui_strcmp(boundaries->nature[izone], "symmetry")) {
2568 inature = CS_SYMMETRY;
2569 }
2570 else if (cs_gui_strcmp(boundaries->nature[izone], "free_inlet_outlet")) {
2571 inature = CS_FREE_INLET;
2572 }
2573 else if (cs_gui_strcmp(boundaries->nature[izone], "free_surface")
2574 && (cs_glob_ale != 0)) {
2575 inature = CS_FREE_SURFACE;
2576 }
2577 else if (cs_gui_strcmp(boundaries->nature[izone], "undefined")) {
2578 inature = CS_INDEF;
2579 }
2580 else if (cs_gui_strcmp(boundaries->nature[izone], "groundwater")) {
2581 inature = CS_INDEF;
2582 }
2583
2584 if (inature < 0)
2585 bft_error(__FILE__, __LINE__, 0,
2586 _("boundary nature %s is unknown \n"),
2587 boundaries->nature[izone]);
2588
2589 int zone_nbr = boundaries->bc_num[izone];
2590
2591 if (nozppm && zone_nbr > *nozppm)
2592 bft_error(__FILE__, __LINE__, 0,
2593 _("zone's label number %i is greater than %i,"
2594 " the maximum allowed \n"), zone_nbr, *nozppm);
2595
2596 } /* for izone */
2597 }
2598
2599 /*============================================================================
2600 * Public function definitions
2601 *============================================================================*/
2602
2603 /*----------------------------------------------------------------------------*/
2604 /*!
2605 * \brief Define boundary conditions based on setup file.
2606 *
2607 * \param[in, out] bdy boundaries structure to update
2608 * (if NULL, default to cs_glob_domain->boundaries)
2609 */
2610 /*----------------------------------------------------------------------------*/
2611
2612 void
2613 cs_gui_boundary_conditions_define(cs_boundary_t *bdy)
2614 {
2615 if (bdy == NULL)
2616 bdy = cs_glob_domain->boundaries;
2617
2618 cs_tree_node_t *tn_b0 = cs_tree_get_node(cs_glob_tree,
2619 "boundary_conditions");
2620
2621 int izone = 0;
2622
2623 /* Wall function info to filter roughness */
2624 cs_wall_functions_t *wall_fnt = cs_get_glob_wall_functions();
2625
2626 /* Build boundary zone definitions */
2627
2628 for (cs_tree_node_t *tn = cs_tree_get_node(tn_b0, "boundary");
2629 tn != NULL;
2630 tn = cs_tree_node_get_next_of_name(tn), izone++) {
2631
2632 /* nature, label and description of the ith boundary zone;
2633 zones are shifted by 1, as default zone 0 is defined
2634 first, before GUI-based definitions (and non-GUI-based
2635 user definitions come last). */
2636
2637 const char *label = cs_tree_node_get_tag(tn, "label");
2638
2639 int bc_num = izone+1;
2640
2641 const int *vi = cs_tree_node_get_child_values_int(tn, "name");
2642 if (vi != NULL)
2643 bc_num = vi[0];
2644
2645 /* label of the ith initialization zone */
2646
2647 const cs_zone_t *z = cs_boundary_zone_by_id(izone + 1);
2648
2649 if (strcmp(label, z->name) != 0)
2650 bft_error(__FILE__, __LINE__, 0,
2651 _("Mismatch between GUI-defined zone %d (%s)\n"
2652 "and boundary condition %d (%s), number %d."),
2653 z->id, z->name, izone+1, label, bc_num);
2654
2655 /* Note: boundary nature is determined again later for most cases,
2656 but at least symmetry is only defined here, and ALE can define
2657 "wall" sections instead of the appropriate type to handle
2658 fixed sections, so also predefine the boundary nature here */
2659
2660 /* Now loop on boundary condition definitions proper */
2661
2662 cs_tree_node_t *tn_bc = NULL;
2663
2664 cs_tree_node_t *tn_b1 = (tn_b0 != NULL) ? tn_b0->children : tn_b0;
2665 for (tn_bc = tn_b1; tn_bc != NULL; tn_bc = tn_bc->next) {
2666
2667 if (cs_gui_strcmp(tn_bc->name, "boundary")) /* handled in parent loop */
2668 continue;
2669
2670 const char *c_label = cs_tree_node_get_tag(tn_bc, "label");
2671 if (c_label != NULL) {
2672 /* Search for label matching boundary */
2673 if (strcmp(c_label, label) == 0)
2674 break;
2675 }
2676
2677 }
2678
2679 if (tn_bc == NULL)
2680 continue;
2681
2682 z = cs_boundary_zone_by_name_try(label);
2683
2684 if (z == NULL) /* may occur when "dead" leaves or present in tree */
2685 continue;
2686
2687 const char *nature = tn_bc->name;
2688
2689 cs_boundary_type_t bc_type = 0;
2690
2691 izone = z->id - 1;
2692 assert(izone >= 0);
2693
2694 if (cs_gui_strcmp(nature, "inlet")) {
2695
2696 bc_type |= CS_BOUNDARY_INLET;
2697
2698 cs_tree_node_t *tn_vp
2699 = cs_tree_node_get_child(tn_bc, "velocity_pressure");
2700
2701 if (cs_glob_physical_model_flag[CS_GROUNDWATER] < 0)
2702 bc_type |= CS_BOUNDARY_IMPOSED_VEL;
2703
2704 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
2705
2706 cs_tree_node_t *tnc = cs_tree_get_node(tn_vp, "compressible_type");
2707 const char *choice = cs_gui_node_get_tag(tnc, "choice");
2708
2709 if (cs_gui_strcmp(choice, "imposed_inlet"))
2710 bc_type |= CS_BOUNDARY_INLET_QH;
2711 else if (cs_gui_strcmp(choice, "subsonic_inlet_PH"))
2712 bc_type |= CS_BOUNDARY_INLET_SUBSONIC_PH;
2713
2714 }
2715
2716 }
2717 else if (cs_gui_strcmp(nature, "wall")) {
2718
2719 bc_type |= CS_BOUNDARY_WALL;
2720
2721 /* sliding wall: velocity */
2722
2723 cs_tree_node_t *tn_vp
2724 = cs_tree_node_get_child(tn_bc, "velocity_pressure");
2725
2726 if (tn_vp != NULL) {
2727 const char *choice = cs_gui_node_get_tag(tn_vp, "choice");
2728 if (cs_gui_strcmp(choice, "on")) {
2729 bc_type |= CS_BOUNDARY_SLIDING_WALL;
2730 _sliding_wall(tn_vp, label);
2731 }
2732
2733 /* check for roughness */
2734 if ( wall_fnt->iwallf != CS_WALL_F_DISABLED
2735 && wall_fnt->iwallf != CS_WALL_F_1SCALE_POWER
2736 && wall_fnt->iwallf != CS_WALL_F_SCALABLE_2SCALES_LOG
2737 && wall_fnt->iwallf != CS_WALL_F_2SCALES_CONTINUOUS) {
2738 cs_real_t roughness = -1.;
2739 cs_gui_node_get_child_real(tn_vp, "roughness", &roughness);
2740 if (roughness > 0) {
2741 bc_type |= CS_BOUNDARY_ROUGH_WALL;
2742 /* Create roughness field if needed */
2743 cs_field_find_or_create("boundary_roughness",
2744 CS_FIELD_INTENSIVE + CS_FIELD_PROPERTY,
2745 CS_MESH_LOCATION_BOUNDARY_FACES,
2746 1, /* dim */
2747 false); /* has previous */
2748 }
2749 }
2750
2751 }
2752 }
2753
2754 else if (cs_gui_strcmp(nature, "outlet")) {
2755
2756 bc_type |= CS_BOUNDARY_OUTLET;
2757
2758 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
2759
2760 cs_tree_node_t *tnc = cs_tree_get_node(tn_bc, "compressible_type");
2761 const char *choice = cs_gui_node_get_tag(tnc, "choice");
2762
2763 if (cs_gui_strcmp(choice, "supersonic_outlet"))
2764 bc_type |= CS_BOUNDARY_SUPERSONIC;
2765 else if (cs_gui_strcmp(choice, "subsonic_outlet"))
2766 bc_type |= CS_BOUNDARY_SUBSONIC;
2767
2768 }
2769
2770 }
2771
2772 else if (cs_gui_strcmp(nature, "free_inlet_outlet")) {
2773 bc_type = bc_type | CS_BOUNDARY_INLET | CS_BOUNDARY_OUTLET;
2774 }
2775
2776 else if (cs_gui_strcmp(nature, "imposed_p_outlet")) {
2777 bc_type = bc_type | CS_BOUNDARY_OUTLET;
2778 bc_type = bc_type | CS_BOUNDARY_IMPOSED_P;
2779 }
2780
2781 else if (!cs_gui_strcmp(nature, "symmetry")) {
2782 bc_type = bc_type | CS_BOUNDARY_SYMMETRY;
2783 }
2784
2785 cs_boundary_add(bdy, bc_type, z->name);
2786 }
2787
2788 /* Definition of the boundaries structure
2789 and some equation parameters */
2790
2791 if (boundaries == NULL)
2792 _init_boundaries();
2793 }
2794
2795 /*----------------------------------------------------------------------------*/
2796 /*!
2797 * \brief Free GUI boundary condition structures.
2798 */
2799 /*----------------------------------------------------------------------------*/
2800
2801 void
2802 cs_gui_boundary_conditions_free_memory(void)
2803 {
2804 int izone;
2805 int n_zones;
2806 int icharb;
2807
2808 /* clean memory for global private structure boundaries */
2809
2810 if (boundaries != NULL) {
2811
2812 n_zones = boundaries->n_zones;
2813
2814 for (int f_id = 0; f_id < boundaries->n_fields; f_id++) {
2815 const cs_field_t *f = cs_field_by_id(f_id);
2816 if (f->type & CS_FIELD_VARIABLE) {
2817 if (boundaries->type_code != NULL)
2818 BFT_FREE(boundaries->type_code[f->id]);
2819 if (boundaries->values != NULL)
2820 BFT_FREE(boundaries->values[f->id]);
2821 if (boundaries->scalar_e != NULL)
2822 BFT_FREE(boundaries->scalar_e[f->id]);
2823 }
2824 }
2825
2826 if ( cs_glob_physical_model_flag[CS_COMBUSTION_PCLC] > -1
2827 || cs_glob_physical_model_flag[CS_COMBUSTION_COAL] > -1) {
2828 const int n_coals = cs_glob_combustion_model->coal.n_coals;
2829 for (izone = 0; izone < n_zones; izone++) {
2830 BFT_FREE(boundaries->qimpcp[izone]);
2831 BFT_FREE(boundaries->timpcp[izone]);
2832 for (icharb = 0; icharb < n_coals; icharb++)
2833 BFT_FREE(boundaries->distch[izone][icharb]);
2834 BFT_FREE(boundaries->distch[izone]);
2835 }
2836 BFT_FREE(boundaries->ientat);
2837 BFT_FREE(boundaries->ientcp);
2838 BFT_FREE(boundaries->inmoxy);
2839 BFT_FREE(boundaries->timpat);
2840 BFT_FREE(boundaries->qimpcp);
2841 BFT_FREE(boundaries->timpcp);
2842 BFT_FREE(boundaries->distch);
2843 }
2844
2845 /* Gas combustion */
2846 BFT_FREE(boundaries->ientfu);
2847 BFT_FREE(boundaries->ientox);
2848 BFT_FREE(boundaries->ientgb);
2849 BFT_FREE(boundaries->ientgf);
2850 BFT_FREE(boundaries->tkent);
2851 BFT_FREE(boundaries->fment);
2852
2853 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
2854 BFT_FREE(boundaries->itype);
2855 BFT_FREE(boundaries->prein);
2856 BFT_FREE(boundaries->rhoin);
2857 BFT_FREE(boundaries->tempin);
2858 }
2859 if (cs_glob_physical_model_flag[CS_GROUNDWATER] > -1) {
2860 BFT_FREE(boundaries->groundwat_e);
2861 }
2862 if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] > -1)
2863 BFT_FREE(boundaries->meteo);
2864
2865 BFT_FREE(boundaries->label);
2866 BFT_FREE(boundaries->nature);
2867 BFT_FREE(boundaries->bc_num);
2868
2869 BFT_FREE(boundaries->iqimp);
2870 BFT_FREE(boundaries->icalke);
2871 BFT_FREE(boundaries->qimp);
2872 BFT_FREE(boundaries->dh);
2873 BFT_FREE(boundaries->xintur);
2874 BFT_FREE(boundaries->type_code);
2875 BFT_FREE(boundaries->values);
2876 BFT_FREE(boundaries->rough);
2877 BFT_FREE(boundaries->norm);
2878 BFT_FREE(boundaries->dir);
2879 BFT_FREE(boundaries->t_to_h);
2880 BFT_FREE(boundaries->velocity_e);
2881 BFT_FREE(boundaries->direction_e);
2882 BFT_FREE(boundaries->scalar_e);
2883 BFT_FREE(boundaries->head_loss_e);
2884
2885 BFT_FREE(boundaries);
2886 }
2887 }
2888
2889 /*----------------------------------------------------------------------------*/
2890
2891 END_C_DECLS
2892