1 /*============================================================================
2 * Management of the GUI parameters file: main parameters
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_all_to_all.h"
52 #include "cs_array.h"
53 #include "cs_base.h"
54 #include "cs_boundary.h"
55 #include "cs_boundary_zone.h"
56 #include "cs_equation.h"
57 #include "cs_equation_param.h"
58 #include "cs_ext_neighborhood.h"
59 #include "cs_field.h"
60 #include "cs_field_pointer.h"
61 #include "cs_file.h"
62 #include "cs_log.h"
63 #include "cs_gui_util.h"
64 #include "cs_gui_boundary_conditions.h"
65 #include "cs_gui_specific_physics.h"
66 #include "cs_gui_mobile_mesh.h"
67 #include "cs_geom.h"
68 #include "cs_gwf_physical_properties.h"
69 #include "cs_internal_coupling.h"
70 #include "cs_math.h"
71 #include "cs_mesh.h"
72 #include "cs_mesh_quantities.h"
73 #include "cs_mesh_location.h"
74 #include "cs_multigrid.h"
75 #include "cs_order.h"
76 #include "cs_parall.h"
77 #include "cs_porous_model.h"
78 #include "cs_parameters.h"
79 #include "cs_param_sles.h"
80 #include "cs_partition.h"
81 #include "cs_physical_model.h"
82 #include "cs_prototypes.h"
83 #include "cs_rotation.h"
84 #include "cs_selector.h"
85 #include "cs_timer.h"
86 #include "cs_time_moment.h"
87 #include "cs_thermal_model.h"
88 #include "cs_physical_properties.h"
89 #include "cs_time_step.h"
90 #include "cs_turbomachinery.h"
91 #include "cs_sles.h"
92 #include "cs_sles_it.h"
93 #include "cs_turbulence_model.h"
94 #include "cs_wall_functions.h"
95 #include "cs_physical_constants.h"
96 #include "cs_balance_by_zone.h"
97 #include "cs_fan.h"
98 #include "cs_velocity_pressure.h"
99 #include "cs_vof.h"
100 #include "cs_volume_zone.h"
101
102 /*----------------------------------------------------------------------------
103 * Header for the current file
104 *----------------------------------------------------------------------------*/
105
106 #include "cs_gui.h"
107
108 /*----------------------------------------------------------------------------*/
109
110 BEGIN_C_DECLS
111
112 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
113
114 /*=============================================================================
115 * Local Macro Definitions
116 *============================================================================*/
117
118 /* debugging switch */
119 #define _XML_DEBUG_ 0
120
121 /*============================================================================
122 * Local Structure Definitions
123 *============================================================================*/
124
125 /*============================================================================
126 * External global variables
127 *============================================================================*/
128
129 /*============================================================================
130 * Private global variables
131 *============================================================================*/
132
133 /*============================================================================
134 * Static global variables
135 *============================================================================*/
136
137 /*============================================================================
138 * Static local variables
139 *============================================================================*/
140
141 /*============================================================================
142 * Private function definitions
143 *============================================================================*/
144
145 /*----------------------------------------------------------------------------
146 * Remove a node inside a tree.
147 *
148 * The node is removed from the tree but not destroyed, so it may be inserted
149 * in another position.
150 *
151 * \param[in, out] tn tree node to remove
152 *----------------------------------------------------------------------------*/
153
154 static void
_tree_node_remove(cs_tree_node_t * tn)155 _tree_node_remove(cs_tree_node_t *tn)
156 {
157 if (tn->prev != NULL)
158 tn->prev->next = tn->next;
159 if (tn->next != NULL)
160 tn->next->prev = tn->prev;
161
162 if (tn->parent != NULL) {
163 if (tn->parent->children == tn)
164 tn->parent->children = tn->next;
165 }
166
167 tn->prev = NULL;
168 tn->next = NULL;
169 }
170
171 /*-----------------------------------------------------------------------------
172 * Modify double numerical parameters.
173 *
174 * parameters:
175 * param <-- label of the numerical parameter
176 * value <-- value of the numerical parameter
177 *----------------------------------------------------------------------------*/
178
179 static void
_numerical_double_parameters(const char * param,double * value)180 _numerical_double_parameters(const char *param,
181 double *value)
182 {
183 cs_tree_node_t *tn = NULL;
184 char *path0 = NULL;
185 BFT_MALLOC(path0, strlen("numerical_parameters/") + strlen(param) + 1, char);
186 strcpy(path0, "numerical_parameters/");
187 strcat(path0, param);
188 tn = cs_tree_get_node(cs_glob_tree, path0);
189 BFT_FREE(path0);
190
191 /* value not changed if path not found */
192 cs_gui_node_get_real(tn, value);
193 }
194
195 /*----------------------------------------------------------------------------
196 * Return the value of the choice attribute for material, method, ...
197 *
198 * parameters:
199 * name <-- name of the property
200 *----------------------------------------------------------------------------*/
201
202 static const char*
_thermal_table_choice(const char * name)203 _thermal_table_choice(const char *name)
204 {
205 cs_tree_node_t *tn
206 = cs_tree_get_node(cs_glob_tree, "physical_properties/fluid_properties");
207 tn = cs_tree_node_get_child(tn, name);
208
209 const char *choice = cs_tree_node_get_tag(tn, "choice");
210
211 return choice;
212 }
213
214 /*----------------------------------------------------------------------------
215 * Return the value of the reference table for a given method, ...
216 *
217 * parameters:
218 * name <-- name of the option
219 *----------------------------------------------------------------------------*/
220
221 static const char*
_thermal_table_option(const char * name)222 _thermal_table_option(const char *name)
223 {
224 cs_tree_node_t *tn
225 = cs_tree_get_node(cs_glob_tree,
226 "physical_properties/fluid_properties/method");
227
228 const char *option = cs_tree_node_get_child_value_str(tn, name);
229
230 return option;
231 }
232
233 /*----------------------------------------------------------------------------
234 * Return the tree node associated to a given property and zone.
235 *----------------------------------------------------------------------------*/
236
237 static cs_tree_node_t *
_get_property_node(const char * property_name,const char * zone_name)238 _get_property_node(const char *property_name,
239 const char *zone_name)
240 {
241 cs_tree_node_t *tn = cs_tree_find_node(cs_glob_tree, "property");
242 while (tn != NULL) {
243 const char *name = cs_tree_node_get_child_value_str(tn, "name");
244 if (cs_gui_strcmp(name, property_name))
245 break;
246 else
247 tn = cs_tree_find_node_next(cs_glob_tree, tn, "property");
248 }
249
250 /* FIXME: Currently zone definitions (other than all_cells) are
251 * handled using sub-nodes. This will be changed for v7.1
252 */
253 if (zone_name != NULL) {
254 if (strcmp(zone_name, "all_cells")) {
255 /* Get zone_id from xml */
256 cs_tree_node_t *id_n =
257 cs_tree_get_node(cs_glob_tree,
258 "solution_domain/volumic_conditions/zone");
259 const char *id_s = NULL;
260 while (id_n != NULL) {
261 const char *zname = cs_tree_node_get_tag(id_n, "label");
262 if (cs_gui_strcmp(zname, zone_name)) {
263 id_s = cs_tree_node_get_tag(id_n, "id");
264 break;
265 }
266 id_n = cs_tree_node_get_next_of_name(id_n);
267 }
268 if (id_s != NULL) {
269 tn = cs_tree_node_get_child(tn, "zone");
270 tn = cs_tree_node_get_sibling_with_tag(tn, "zone_id", id_s);
271 }
272 }
273 }
274
275 return tn;
276 }
277
278 /*----------------------------------------------------------------------------
279 * Return the value of the choice attribute from a property name.
280 *
281 * parameters:
282 * property_name <-- name of the property
283 * zone_name <-- name of zone. If NULL we use
284 * default definition
285 *----------------------------------------------------------------------------*/
286
287 static const char*
_properties_choice(const char * property_name,const char * zone_name)288 _properties_choice(const char *property_name,
289 const char *zone_name)
290 {
291 const char *choice = NULL;
292
293 cs_tree_node_t *node = _get_property_node(property_name, zone_name);
294
295 choice = cs_tree_node_get_child_value_str(node, "choice");
296
297 return choice;
298 }
299
300 /*----------------------------------------------------------------------------
301 * Return the formula associated to a given property and zone.
302 *----------------------------------------------------------------------------*/
303
304 static const char*
_property_formula(const char * property_name,const char * zone_name)305 _property_formula(const char *property_name,
306 const char *zone_name)
307 {
308 const char *law = NULL;
309
310 cs_tree_node_t *node = _get_property_node(property_name, zone_name);
311
312 node = cs_tree_node_get_child(node, "formula");
313
314 law = cs_tree_node_get_value_str(node);
315
316 return law;
317 }
318
319 /*----------------------------------------------------------------------------
320 * Return 0 if default value is needed
321 *
322 * parameters:
323 * name <-- name of the property
324 *----------------------------------------------------------------------------*/
325
326 static int
_thermal_table_needed(const char * name)327 _thermal_table_needed(const char *name)
328 {
329 int choice = 0;
330
331 const char *prop_choice = _properties_choice(name, NULL);
332 if (cs_gui_strcmp(prop_choice, "thermal_law"))
333 choice = 1;
334
335 return choice;
336 }
337
338 /*-----------------------------------------------------------------------------
339 * Compute a physical property based on a thermal law.
340 *----------------------------------------------------------------------------*/
341
342 static void
_physical_property_thermal_law(cs_field_t * c_prop,const cs_zone_t * z,cs_phys_prop_type_t property)343 _physical_property_thermal_law(cs_field_t *c_prop,
344 const cs_zone_t *z,
345 cs_phys_prop_type_t property)
346 {
347 /* For incompressible flows, the thermodynamic pressure is constant over
348 * time and is the reference pressure. */
349
350 cs_lnum_t thermodynamic_pressure_stride = 0;
351 cs_lnum_t thermal_f_val_stride = 1;
352 cs_real_t _p0 = cs_glob_fluid_properties->p0;
353 cs_real_t _t0 = cs_glob_fluid_properties->t0;
354
355 const cs_real_t *thermodynamic_pressure = &_p0;
356 const cs_real_t *_thermal_f_val = NULL;
357
358 if (CS_F_(t) != NULL) {
359 if (CS_F_(t)->type & CS_FIELD_VARIABLE)
360 _thermal_f_val = CS_F_(t)->val;
361 }
362 else if (CS_F_(h) != NULL) {
363 if (CS_F_(h)->type & CS_FIELD_VARIABLE)
364 _thermal_f_val = CS_F_(h)->val;
365 }
366 else if (CS_F_(e_tot) != NULL) {
367 if (CS_F_(h)->type & CS_FIELD_VARIABLE) {
368 _thermal_f_val = CS_F_(e_tot)->val;
369 thermodynamic_pressure = CS_F_(p)->val;
370 thermodynamic_pressure_stride = 1;
371 }
372 }
373 else {
374 thermal_f_val_stride = 0;
375 _thermal_f_val = &_t0;
376 }
377
378 cs_phys_prop_compute(property,
379 z->n_elts,
380 thermodynamic_pressure_stride,
381 thermal_f_val_stride,
382 thermodynamic_pressure,
383 _thermal_f_val,
384 c_prop->val);
385 }
386
387 /*-----------------------------------------------------------------------------
388 * use MEI for thermal diffusivity
389 *
390 * This is a special case of _physical_property_, where the thermal
391 * diffusivity computation is based on that of the thermal conductivity.
392 *----------------------------------------------------------------------------*/
393
394 static void
_physical_property_th_diffusivity(cs_field_t * c_prop,const cs_zone_t * z)395 _physical_property_th_diffusivity(cs_field_t *c_prop,
396 const cs_zone_t *z)
397 {
398 const char prop_name[] = "thermal_conductivity";
399
400 const cs_lnum_t n_cells = z->n_elts;
401 const cs_lnum_t *cell_ids = z->elt_ids;
402
403 int user_law = 0;
404 const char *prop_choice = _properties_choice(prop_name, NULL);
405
406 if (cs_gui_strcmp(prop_choice, "user_law"))
407 user_law = 1;
408 else if (z->id > 1 && z->type & CS_VOLUME_ZONE_PHYSICAL_PROPERTIES) {
409 const char *law = _property_formula(prop_name, z->name);
410 if (law != NULL)
411 user_law = 1;
412 }
413
414 if (user_law) {
415
416 const char *law = _property_formula(prop_name, z->name);
417
418 if (law != NULL) {
419 /* Pass "field overlay": shallow copy of field with modified name
420 so as to be handled by the MEG volume function. */
421 cs_field_t _c_prop = *c_prop;
422 _c_prop.name = prop_name;
423 cs_field_t *fmeg[1] = {&_c_prop};
424 cs_meg_volume_function(z, fmeg);
425 }
426
427 }
428 else if (cs_gui_strcmp(prop_choice, "thermal_law") &&
429 cs_gui_strcmp(z->name, "all_cells")) {
430
431 _physical_property_thermal_law(c_prop,
432 z,
433 CS_PHYS_PROP_THERMAL_CONDUCTIVITY);
434
435 }
436
437 /* If property is globaly variable but defined as a constant
438 for this zone, set from constant */
439
440 else if (z->id <= 1 && cs_gui_strcmp(prop_choice, "constant")) {
441
442 const cs_fluid_properties_t *phys_pp = cs_glob_fluid_properties;
443 const cs_real_t p_val = phys_pp->lambda0;
444
445 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
446 cs_lnum_t c_id = cell_ids[e_id];
447 c_prop->val[c_id] = p_val;
448 }
449
450 }
451
452 /* Finalize special case for conduction to diffusion conversion */
453
454 if (CS_F_(cp) == NULL) {
455 const cs_real_t cp0 = cs_glob_fluid_properties->cp0;
456 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
457 cs_lnum_t c_id = cell_ids[e_id];
458 c_prop->val[c_id] /= cp0;
459 }
460 }
461 else {
462 const cs_real_t *cp = CS_F_(cp)->val;
463 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
464 cs_lnum_t c_id = cell_ids[e_id];
465 c_prop->val[c_id] /= cp[c_id];
466 }
467 }
468 }
469
470 /*-----------------------------------------------------------------------------
471 * use MEI for physical property
472 *----------------------------------------------------------------------------*/
473
474 static void
_physical_property(cs_field_t * c_prop,const cs_zone_t * z)475 _physical_property(cs_field_t *c_prop,
476 const cs_zone_t *z)
477 {
478 int user_law = 0;
479 const char *prop_choice = _properties_choice(c_prop->name, NULL);
480
481 /* Special case: "thermal_conductivity" is defined by GUI, but
482 in case of Enthalpy thermal model, the matching field is
483 "thermal_diffusivity". */
484
485 if (prop_choice == NULL) {
486 if ( cs_glob_thermal_model->itherm == CS_THERMAL_MODEL_ENTHALPY
487 && cs_gui_strcmp(c_prop->name, "thermal_diffusivity")) {
488 _physical_property_th_diffusivity(c_prop, z);
489 return;
490 }
491 }
492
493 if (cs_gui_strcmp(prop_choice, "user_law"))
494 user_law = 1;
495 else if (z->id > 1 && z->type & CS_VOLUME_ZONE_PHYSICAL_PROPERTIES) {
496 const char *law = _property_formula(c_prop->name, z->name);
497 if (law != NULL)
498 user_law = 1;
499 }
500
501 if (user_law) {
502
503 const char *law = _property_formula(c_prop->name, z->name);
504
505 if (law != NULL) {
506 cs_field_t *fmeg[1] = {c_prop};
507 cs_meg_volume_function(z, fmeg);
508 }
509
510 }
511 else if (cs_gui_strcmp(prop_choice, "thermal_law") &&
512 cs_gui_strcmp(z->name, "all_cells")) {
513 cs_phys_prop_type_t property = -1;
514
515 if (cs_gui_strcmp(c_prop->name, "density"))
516 property = CS_PHYS_PROP_DENSITY;
517
518 else if (cs_gui_strcmp(c_prop->name, "molecular_viscosity"))
519 property = CS_PHYS_PROP_DYNAMIC_VISCOSITY;
520
521 else if (cs_gui_strcmp(c_prop->name, "specific_heat"))
522 property = CS_PHYS_PROP_ISOBARIC_HEAT_CAPACITY;
523
524 else if (cs_gui_strcmp(c_prop->name, "thermal_conductivity"))
525 property = CS_PHYS_PROP_THERMAL_CONDUCTIVITY;
526
527 else
528 bft_error(__FILE__, __LINE__, 0,
529 _("Error: can not evaluate property: %s using a thermal law\n"),
530 c_prop->name);
531
532 _physical_property_thermal_law(c_prop, z, property);
533 }
534
535 /* If property is globaly variable but defined as a constant
536 for this zone, set from constant */
537
538 else if (z->id <= 1 && cs_gui_strcmp(prop_choice, "constant")) {
539 cs_real_t p_val = -1;
540 const cs_fluid_properties_t *phys_pp = cs_glob_fluid_properties;
541
542 if (cs_gui_strcmp(c_prop->name, "density"))
543 p_val = phys_pp->ro0;
544 else if (cs_gui_strcmp(c_prop->name, "molecular_viscosity"))
545 p_val = phys_pp->viscl0;
546 else if (cs_gui_strcmp(c_prop->name, "specific_heat"))
547 p_val = phys_pp->cp0;
548 else if (cs_gui_strcmp(c_prop->name, "thermal_conductivity")) {
549 p_val = phys_pp->lambda0;
550 }
551
552 const cs_lnum_t n_cells = z->n_elts;
553 const cs_lnum_t *cell_ids = z->elt_ids;
554 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
555 cs_lnum_t c_id = cell_ids[e_id];
556 c_prop->val[c_id] = p_val;
557 }
558 }
559 }
560
561 /*-----------------------------------------------------------------------------
562 * Return the value of choice for user scalar's property
563 *
564 * parameters:
565 * f_name <-- field name
566 * choice <-> choice for property
567 *----------------------------------------------------------------------------*/
568
569 static int
_scalar_properties_choice(const char * f_name,int * choice)570 _scalar_properties_choice(const char *f_name,
571 int *choice)
572 {
573 const char *buff = NULL;
574 int ichoice;
575
576 cs_tree_node_t *tn;
577
578 for (tn = cs_tree_get_node(cs_glob_tree, "additional_scalars/variable");
579 tn != NULL;) {
580 if (!cs_gui_strcmp(f_name, cs_tree_node_get_tag(tn, "name")))
581 tn = cs_tree_node_get_next_of_name(tn);
582 else
583 break;
584 }
585
586 tn = cs_tree_get_node(tn, "property/choice");
587 buff = cs_tree_node_get_value_str(tn);
588
589 if (buff == NULL) {
590 ichoice = 0;
591
592 } else {
593 ichoice = 1;
594
595 if (cs_gui_strcmp(buff, "user_law"))
596 *choice = 1;
597 else if (cs_gui_strcmp(buff, "constant"))
598 *choice = 0;
599 else
600 bft_error(__FILE__, __LINE__, 0, _("Invalid node in function %s\n"),
601 __func__);
602 }
603
604 return ichoice;
605 }
606
607 /*-----------------------------------------------------------------------------
608 * Return value of diffusion coefficient for user scalars
609 * return 1 if value exists
610 * return 0 if not
611 *
612 * parameters:
613 * f <-- pointer to field
614 * value <-- value of diffusion coefficient
615 *----------------------------------------------------------------------------*/
616
617 static void
_scalar_diffusion_value(const cs_field_t * f,cs_real_t * value)618 _scalar_diffusion_value(const cs_field_t *f,
619 cs_real_t *value)
620 {
621 cs_tree_node_t *tn
622 = cs_tree_get_node(cs_glob_tree, "additional_scalars/variable");
623
624 while (tn != NULL) {
625 const char *name = cs_tree_node_get_child_value_str(tn, "name");
626 if (cs_gui_strcmp(name, f->name))
627 break;
628 else
629 tn = cs_tree_find_node_next(cs_glob_tree, tn, "variable");
630 }
631
632 tn = cs_tree_get_node(tn, "property/initial_value");
633
634 cs_gui_node_get_real(tn, value);
635 }
636
637 /*-----------------------------------------------------------------------------
638 * Find the node associated with a given variable.
639 *
640 * parameters:
641 * variable_name <-- name of variable
642 *
643 * return:
644 * pointer to associated tree node (NULL if not present)
645 *----------------------------------------------------------------------------*/
646
647 static cs_tree_node_t *
_find_node_variable(const char * variable_name)648 _find_node_variable(const char *variable_name)
649 {
650 cs_tree_node_t *tn = cs_tree_find_node(cs_glob_tree, "variable");
651 while (tn != NULL) {
652 const char *name = cs_tree_node_get_child_value_str(tn, "name");
653 if (cs_gui_strcmp(name, variable_name))
654 break;
655 else
656 tn = cs_tree_find_node_next(cs_glob_tree, tn, "variable");
657 }
658
659 return tn;
660 }
661
662 /*-----------------------------------------------------------------------------
663 * Return value of turbulent flux model
664 *
665 * parameters:
666 * tn_v <-- tree node associated with variable
667 * value --> value of turbulent flux model
668 *----------------------------------------------------------------------------*/
669
670 static void
_variable_turbulent_flux_model(cs_tree_node_t * tn_v,int * value)671 _variable_turbulent_flux_model(cs_tree_node_t *tn_v,
672 int *value)
673 {
674 const char *result
675 = cs_tree_node_get_child_value_str(tn_v, "turbulent_flux_model");
676
677 if (cs_gui_strcmp(result, "SGDH"))
678 *value = 0;
679 else if (cs_gui_strcmp(result, "GGDH"))
680 *value = 10;
681 else if (cs_gui_strcmp(result, "EB-GGDH"))
682 *value = 11;
683 else if (cs_gui_strcmp(result, "AFM"))
684 *value = 20;
685 else if (cs_gui_strcmp(result, "EB-AFM"))
686 *value = 21;
687 else if (cs_gui_strcmp(result, "DFM"))
688 *value = 30;
689 else if (cs_gui_strcmp(result, "EB-DFM"))
690 *value = 31;
691 else
692 *value = 0; /* assign default */
693 }
694
695 /*----------------------------------------------------------------------------
696 * Get the attribute value from the scheme order
697 *
698 * parameters:
699 * tn_v <-- node assocaited with variable
700 * keyword --> value of attribute node
701 *----------------------------------------------------------------------------*/
702
703 static void
_order_scheme_value(cs_tree_node_t * tn_v,int * keyword)704 _order_scheme_value(cs_tree_node_t *tn_v,
705 int *keyword)
706 {
707 cs_tree_node_t *tn = cs_tree_node_get_child(tn_v, "order_scheme");
708 const char *choice = cs_tree_node_get_child_value_str(tn, "choice");
709
710 if (cs_gui_strcmp(choice, "centered"))
711 *keyword = 1;
712 else if (cs_gui_strcmp(choice, "solu"))
713 *keyword = 0;
714 }
715
716 /*----------------------------------------------------------------------------
717 * Get the attribute value from the cs_tree query.
718 *
719 * parameters:
720 * node <--
721 * child <-- child markup
722 * keyword --> value of attribute node
723 *----------------------------------------------------------------------------*/
724
725 static void
_slope_test_value(cs_tree_node_t * tn,int * keyword)726 _slope_test_value(cs_tree_node_t *tn,
727 int *keyword)
728 {
729 int result = -999;
730 cs_gui_node_get_child_status_int(tn, "slope_test", &result);
731
732 if (result == 1)
733 *keyword = 0;
734 if (result == 0)
735 *keyword = 1;
736 }
737
738 /*----------------------------------------------------------------------------
739 * Return the attribute choice associated to a child markup from a variable.
740 *
741 * parameters:
742 * tn_v <-- tree node associated with the choice
743 * child <-- child markup
744 *----------------------------------------------------------------------------*/
745
746 static const char *
_variable_choice(cs_tree_node_t * tn_v,const char * child)747 _variable_choice(cs_tree_node_t *tn_v,
748 const char *child)
749 {
750 cs_tree_node_t *tn = cs_tree_get_node(tn_v, child);
751
752 const char *choice = cs_tree_node_get_child_value_str(tn, "choice");
753
754 return choice;
755 }
756
757 /*-----------------------------------------------------------------------------
758 * Modify integer numerical parameters.
759 *
760 * parameters:
761 * param <-- label of the numerical parameter
762 * keyword <-- value of the numerical parameter
763 *----------------------------------------------------------------------------*/
764
765 static void
_numerical_int_parameters(const char * param,int * keyword)766 _numerical_int_parameters(const char *param,
767 int *keyword)
768 {
769 int choice = *keyword;
770 int result = *keyword;
771 cs_tree_node_t *tn = cs_tree_get_node(cs_glob_tree, "numerical_parameters");
772
773 if (cs_gui_strcmp(param, "gradient_reconstruction")) {
774
775 tn = cs_tree_get_node(tn, param);
776 tn = cs_tree_get_node(tn, "choice");
777 cs_gui_node_get_int(tn, &choice);
778 *keyword = choice;
779
780 } else if (cs_gui_strcmp(param, "piso_sweep_number")) {
781
782 tn = cs_tree_get_node(tn, "velocity_pressure_algo");
783 tn = cs_tree_get_node(tn, param);
784 cs_gui_node_get_int(tn, &result);
785 *keyword = result;
786
787 } else {
788
789 tn = cs_tree_get_node(tn, param);
790 cs_gui_node_get_status_int(tn, &result);
791 *keyword = result;
792 }
793 }
794
795 /*-----------------------------------------------------------------------------
796 * Modify gravity parameters.
797 *
798 * parameters:
799 * param <-- gravity parameter (GX, GY, GZ)
800 * keyword <<-- new value of the gravity parameter
801 *----------------------------------------------------------------------------*/
802
803 static void
_gravity_value(const char * param,double * value)804 _gravity_value(const char *param,
805 double *value)
806 {
807 cs_tree_node_t *tn
808 = cs_tree_get_node(cs_glob_tree, "physical_properties/gravity");
809 tn = cs_tree_get_node(tn, param);
810
811 cs_gui_node_get_real(tn, value);
812 }
813
814 /*-----------------------------------------------------------------------------
815 * Modify coriolis source terms parameters.
816 *
817 * parameters:
818 * param <-- coriolis parameter (omegax, omegay, omegaz)
819 * value --> new value of the coriolis parameter
820 *----------------------------------------------------------------------------*/
821
822 static void
_coriolis_value(const char * param,double * value)823 _coriolis_value(const char *param,
824 double *value)
825 {
826 cs_tree_node_t *tn
827 = cs_tree_get_node(cs_glob_tree, "physical_properties/omega");
828 tn = cs_tree_get_node(tn, param);
829
830 cs_gui_node_get_real(tn, value);
831 }
832
833 /*----------------------------------------------------------------------------
834 * Get the value of the choice attribute from a property markup.
835 * Return 1 if the cs_tree request has succeeded, 0 otherwise.
836 *
837 * parameters:
838 * property_name <-- name of the property
839 * choice --> value of the attribute choice
840 *----------------------------------------------------------------------------*/
841
842 static int
_properties_choice_id(const char * property_name,int * choice)843 _properties_choice_id(const char *property_name,
844 int *choice)
845 {
846 const char *buff = NULL;
847 int iok = 0;
848
849 buff = _properties_choice(property_name, NULL);
850 *choice = 0; /* default */
851 if (buff)
852 {
853 iok = 1;
854 if (cs_gui_strcmp(buff, "user_law") ||
855 cs_gui_strcmp(buff, "predefined_law") ||
856 cs_gui_strcmp(buff, "thermal_law") )
857 *choice = 1;
858 else if (cs_gui_strcmp(buff, "constant"))
859 *choice = 0;
860 }
861 else
862 iok = 0;
863 return iok;
864 }
865
866 /*----------------------------------------------------------------------------
867 * Return the length choice for initialize turbulence
868 *----------------------------------------------------------------------------*/
869
870 static const char *
_reference_length_initialization_choice(void)871 _reference_length_initialization_choice(void)
872 {
873 cs_tree_node_t *tn
874 = cs_tree_get_node(cs_glob_tree,
875 "thermophysical_models/turbulence/"
876 "reference_length/choice");
877
878 const char *initialization_choice
879 = cs_tree_node_get_value_str(tn);
880
881 return initialization_choice;
882 }
883
884 /*----------------------------------------------------------------------------
885 * Return the initialization choice of the turbulence variables.
886 *
887 * parameters:
888 * zone_id <-- zone number
889 *----------------------------------------------------------------------------*/
890
891 static const char *
_turbulence_initialization_choice(const char * zone_id)892 _turbulence_initialization_choice(const char* zone_id)
893 {
894 cs_tree_node_t *tn
895 = cs_tree_get_node(cs_glob_tree,
896 "thermophysical_models/turbulence/initialization");
897 tn = cs_tree_node_get_sibling_with_tag(tn, "zone_id", zone_id);
898 tn = cs_tree_get_node(tn, "choice");
899
900 const char *initialization_choice
901 = cs_tree_node_get_value_str(tn);
902
903 return initialization_choice;
904 }
905
906 /*==========================
907 * FOR VOLUMIC ZONES
908 *==========================*/
909
910 /*----------------------------------------------------------------------------
911 * Check if a given zone type attribute is active
912 *
913 * parameters:
914 * z_id <-- zone id
915 * attr <-- attribute checked for
916 *
917 * return:
918 * true if attribute checked for is active, false otherwise
919 *----------------------------------------------------------------------------*/
920
921 static bool
_zone_id_is_type(int z_id,const char * attr)922 _zone_id_is_type(int z_id,
923 const char *attr)
924 {
925 bool retval = false;
926
927 const char *status = NULL;
928 cs_tree_node_t *tn
929 = cs_tree_get_node(cs_glob_tree, "solution_domain/volumic_conditions/zone");
930 for (int i = 1; tn != NULL && i < z_id; i++) {
931 tn = cs_tree_node_get_next_of_name(tn);
932 }
933 tn = cs_tree_get_node(tn, attr);
934 status = cs_tree_node_get_value_str(tn);
935
936 if (status != NULL) {
937 if (cs_gui_strcmp(status, "on"))
938 retval = true;
939 }
940
941 return retval;
942 }
943
944 /*----------------------------------------------------------------------------
945 * Indicate if a given zone type attribute is active
946 *
947 * parameters:
948 * tn <-- tree node associated with the zone variable
949 * type_str <-- string describing the type
950 *----------------------------------------------------------------------------*/
951
952 static bool
_zone_is_type(cs_tree_node_t * tn,const char * type_str)953 _zone_is_type(cs_tree_node_t *tn,
954 const char *type_str)
955 {
956 bool retval = false;
957
958 const char *type_s = cs_tree_node_get_tag(tn, type_str);
959 if (type_s != NULL) {
960 if (strcmp(type_s, "on") == 0)
961 retval = true;
962 }
963
964 return retval;
965 }
966
967 /*----------------------------------------------------------------------------
968 * Add a test attribute to an cs_tree query for a given zone id.
969 *
970 * This should dissapear in the future, when zone names replace ids
971 * in the XML file.
972 *
973 * Note that ids in the zone model are 0-based, while those in the XML
974 * are 1-based. The shift is handled by adding a default zone in
975 * the base model.
976 *
977 * parameters:
978 * tn <-- parent tree node
979 * z_id <-- zone id (O-based)
980 *
981 * return:
982 * sibling node matching attribute
983 *----------------------------------------------------------------------------*/
984
985 static cs_tree_node_t *
_add_zone_id_test_attribute(cs_tree_node_t * tn,int z_id)986 _add_zone_id_test_attribute(cs_tree_node_t *tn,
987 int z_id)
988 {
989 char z_id_str[32];
990 snprintf(z_id_str, 31, "%d", z_id);
991
992 return cs_tree_node_get_sibling_with_tag(tn, "zone_id", z_id_str);
993 }
994
995 /*-----------------------------------------------------------------------------
996 * Change the head losses matrix from the local frame to the global frame.
997 *
998 * parameters:
999 * a_ij <-- change matrix from the local frame to the global frame
1000 * in_ij <-- head losses matrix in the local frame
1001 * out_ij <-- head losses matrix in the global frame
1002 *----------------------------------------------------------------------------*/
1003
1004 static void
_matrix_base_conversion(double a11,double a12,double a13,double a21,double a22,double a23,double a31,double a32,double a33,double in11,double in12,double in13,double in21,double in22,double in23,double in31,double in32,double in33,double * out11,double * out12,double * out13,double * out21,double * out22,double * out23,double * out31,double * out32,double * out33)1005 _matrix_base_conversion(double a11, double a12, double a13,
1006 double a21, double a22, double a23,
1007 double a31, double a32, double a33,
1008 double in11, double in12, double in13,
1009 double in21, double in22, double in23,
1010 double in31, double in32, double in33,
1011 double *out11, double *out12, double *out13,
1012 double *out21, double *out22, double *out23,
1013 double *out31, double *out32, double *out33)
1014 {
1015 int i, j, k;
1016 double tensorP[3][3], tensorA[3][3], tensorB[3][3];
1017 double tensorC[3][3], tensorD[3][3];
1018
1019 tensorA[0][0] = in11;
1020 tensorA[0][1] = in12;
1021 tensorA[0][2] = in13;
1022 tensorA[1][0] = in21;
1023 tensorA[1][1] = in22;
1024 tensorA[1][2] = in23;
1025 tensorA[2][0] = in31;
1026 tensorA[2][1] = in32;
1027 tensorA[2][2] = in33;
1028
1029 tensorP[0][0] = a11;
1030 tensorP[0][1] = a12;
1031 tensorP[0][2] = a13;
1032 tensorP[1][0] = a21;
1033 tensorP[1][1] = a22;
1034 tensorP[1][2] = a23;
1035 tensorP[2][0] = a31;
1036 tensorP[2][1] = a32;
1037 tensorP[2][2] = a33;
1038
1039 for (i = 0; i < 3; i++) {
1040 for (j = 0; j < 3; j++) {
1041 tensorB[i][j] = 0.;
1042 for (k = 0; k < 3; k++)
1043 tensorB[i][j] += tensorP[i][k] * tensorA[k][j];
1044 }
1045 }
1046
1047 /* Inversion of a 3x3 matrix */
1048
1049 tensorC[0][0] = a11;
1050 tensorC[0][1] = a21;
1051 tensorC[0][2] = a31;
1052 tensorC[1][0] = a12;
1053 tensorC[1][1] = a22;
1054 tensorC[1][2] = a32;
1055 tensorC[2][0] = a13;
1056 tensorC[2][1] = a23;
1057 tensorC[2][2] = a33;
1058
1059 for (i = 0; i < 3; i++) {
1060 for (j = 0; j < 3; j++) {
1061 tensorD[i][j] = 0.;
1062 for (k = 0; k < 3; k++)
1063 tensorD[i][j] += tensorB[i][k] * tensorC[k][j];
1064 }
1065 }
1066
1067 *out11 = tensorD[0][0];
1068 *out22 = tensorD[1][1];
1069 *out33 = tensorD[2][2];
1070 *out12 = tensorD[0][1];
1071 *out13 = tensorD[0][2];
1072 *out21 = tensorD[1][0];
1073 *out23 = tensorD[1][2];
1074 *out31 = tensorD[2][0];
1075 *out32 = tensorD[2][1];
1076 }
1077
1078 /*-----------------------------------------------------------------------------
1079 * Return value of coefficient associated to the head losses definition.
1080 *
1081 * parameters:
1082 * tn_hl <-- tree node associated with zone head losses
1083 * c <-- name of the coefficient
1084 *----------------------------------------------------------------------------*/
1085
1086 static double
_c_head_losses(cs_tree_node_t * tn_hl,const char * c)1087 _c_head_losses(cs_tree_node_t *tn_hl,
1088 const char *c)
1089 {
1090 double value = 0.0;
1091
1092 const cs_real_t *v_r = cs_tree_node_get_child_values_real(tn_hl, c);
1093 if (v_r != NULL)
1094 value = v_r[0];
1095
1096 return value;
1097 }
1098
1099 /*-----------------------------------------------------------------------------
1100 * Get turbomachinery model
1101 *
1102 * parameters:
1103 * model_type --> turbomachinery model type
1104 * coupled --> use coupled variant
1105 *----------------------------------------------------------------------------*/
1106
1107 static void
_turbomachinery_model(cs_turbomachinery_model_t * model_type,bool * coupled)1108 _turbomachinery_model(cs_turbomachinery_model_t *model_type,
1109 bool *coupled)
1110 {
1111 *model_type = CS_TURBOMACHINERY_NONE;
1112 *coupled = false;
1113
1114 const char *model = NULL;
1115
1116 cs_tree_node_t *tn
1117 = cs_tree_get_node(cs_glob_tree,
1118 "thermophysical_models/turbomachinery/model");
1119 model = cs_tree_node_get_value_str(tn);
1120
1121 if (cs_gui_strcmp(model, "off"))
1122 *model_type = CS_TURBOMACHINERY_NONE;
1123 else if (cs_gui_strcmp(model, "transient"))
1124 *model_type = CS_TURBOMACHINERY_TRANSIENT;
1125 else if (cs_gui_strcmp(model, "frozen"))
1126 *model_type = CS_TURBOMACHINERY_FROZEN;
1127 else if (cs_gui_strcmp(model, "transient_coupled")) {
1128 *model_type = CS_TURBOMACHINERY_TRANSIENT;
1129 *coupled = true;
1130 }
1131 else if (cs_gui_strcmp(model, "frozen_coupled")) {
1132 *model_type = CS_TURBOMACHINERY_FROZEN;
1133 *coupled = true;
1134 }
1135 }
1136
1137 /*----------------------------------------------------------------------------
1138 * Return the value of the choice attribute for rotor (turbomachinery)
1139 *
1140 * parameters:
1141 * rotor_id <-- id of the rotor
1142 * name <-- name of the property
1143 *----------------------------------------------------------------------------*/
1144
1145 static double
_rotor_option(int rotor_id,const char * name)1146 _rotor_option(int rotor_id,
1147 const char *name)
1148 {
1149 double value = 0.;
1150
1151 cs_tree_node_t *tn
1152 = cs_tree_get_node(cs_glob_tree,
1153 "thermophysical_models/turbomachinery/rotor");
1154 for (int i = 1; tn != NULL && i < rotor_id + 1; i++) {
1155 tn = cs_tree_node_get_next_of_name(tn);
1156 }
1157 tn = cs_tree_node_get_child(tn, "rotation");
1158 tn = cs_tree_node_get_child(tn, name);
1159
1160 cs_gui_node_get_real(tn, &value);
1161
1162 return value;
1163 }
1164
1165 /*-----------------------------------------------------------------------------
1166 * Return the value to a face joining markup for turbomachinery
1167 *
1168 * parameters:
1169 * keyword <-- label of the markup
1170 * number <-- joining number
1171 *----------------------------------------------------------------------------*/
1172
1173 static const char *
_get_rotor_face_joining(const char * keyword,int number)1174 _get_rotor_face_joining(const char *keyword,
1175 int number)
1176 {
1177 cs_tree_node_t *tn
1178 = cs_tree_get_node(cs_glob_tree,
1179 "thermophysical_models/turbomachinery/"
1180 "joining/face_joining");
1181 for (int i = 1; tn != NULL && i < number; i++) {
1182 tn = cs_tree_node_get_next_of_name(tn);
1183 }
1184 tn = cs_tree_get_node(tn, keyword);
1185
1186 const char *value = cs_tree_node_get_value_str(tn);
1187
1188 return value;
1189 }
1190
1191 /*----------------------------------------------------------------------------
1192 * Return volume zone id in tree (1 to n).
1193 *
1194 * Also checks that zone definitions are ordered if id_e > -1
1195 *
1196 * Note that the name tag for zones actually defines an integer (1 to n).
1197 * This tag should be removed in the future to only use the zone label
1198 * (the actual zone name).
1199 *
1200 * parameters:
1201 * tn <-- associated tree node
1202 * id_e <-- expected zone id (0 to n-1)
1203 *
1204 * return
1205 * zone's id in tree (1 to n)
1206 *----------------------------------------------------------------------------*/
1207
1208 static int
_v_zone_t_id(cs_tree_node_t * tn,int id_e)1209 _v_zone_t_id(cs_tree_node_t *tn,
1210 int id_e)
1211 {
1212 int z_t_id = id_e + 1;
1213
1214 const char *id_s = cs_tree_node_get_tag(tn, "id");
1215 if (id_s != NULL) {
1216 z_t_id = atoi(id_s);
1217 if (id_e > -1 && z_t_id != id_e + 1)
1218 bft_printf(_("\n"
1219 " Warning: noncontiguous %s zone ids in XML:\n"
1220 " zone with index %d has id %d.\n"),
1221 tn->name, id_e, z_t_id);
1222 }
1223
1224 return z_t_id;
1225 }
1226
1227 /*----------------------------------------------------------------------------
1228 * Return volume zone tree node with a given id (name)
1229 *
1230 * parameters:
1231 * tn_vc <-- associated volume conditions tree node (parent)
1232 * z_t_id <-- zone id in tree (1 to n)
1233 *
1234 * return
1235 * pointer to associated zone node, or NULL if not found
1236 *----------------------------------------------------------------------------*/
1237
1238 static cs_tree_node_t *
_v_zone_node_by_id(cs_tree_node_t * tn_vc,int z_t_id)1239 _v_zone_node_by_id(cs_tree_node_t *tn_vc,
1240 int z_t_id)
1241 {
1242 cs_tree_node_t *retval = NULL;
1243
1244 for (cs_tree_node_t *tn = cs_tree_node_get_child(tn_vc, "zone");
1245 tn != NULL;
1246 tn = cs_tree_node_get_next_of_name(tn)) {
1247 if (z_t_id == _v_zone_t_id(tn, -1)) {
1248 retval = tn;
1249 break;
1250 }
1251 }
1252 return retval;
1253 }
1254
1255 /*----------------------------------------------------------------------------
1256 * Ensure volume and boundary zones defined through the GUI are ordered
1257 * (modifying tree if necessary)
1258 *----------------------------------------------------------------------------*/
1259
1260 static void
_ensure_zones_order(void)1261 _ensure_zones_order(void)
1262 {
1263 cs_tree_node_t *tn_parent = NULL;
1264
1265 /* Volume zones */
1266 /*------------- */
1267
1268 tn_parent = cs_tree_get_node(cs_glob_tree,
1269 "solution_domain/volumic_conditions");
1270
1271 /* Check if volume zones are defined in increasing id order */
1272
1273 bool need_reorder = false;
1274
1275 int z_id_prev = -1;
1276 int id = 0;
1277
1278 for (cs_tree_node_t *tn = cs_tree_node_get_child(tn_parent, "zone");
1279 tn != NULL;
1280 tn = cs_tree_node_get_next_of_name(tn), id++) {
1281 if (_v_zone_t_id(tn, id) < z_id_prev)
1282 need_reorder = true;
1283 }
1284
1285 const int n_v_zones = id;
1286
1287 if (need_reorder) {
1288
1289 cs_lnum_t *order = NULL, *z_ids = NULL;
1290
1291 /* Build ordering array */
1292
1293 BFT_MALLOC(z_ids, n_v_zones, cs_lnum_t);
1294 BFT_MALLOC(order, n_v_zones, cs_lnum_t);
1295
1296 /* Loop on volume condition zones */
1297
1298 id = 0;
1299 for (cs_tree_node_t *tn = cs_tree_node_get_child(tn_parent, "zone");
1300 tn != NULL;
1301 tn = cs_tree_node_get_next_of_name(tn), id++) {
1302 z_ids[id] = _v_zone_t_id(tn, id);
1303 }
1304
1305 cs_order_lnum_allocated(NULL, z_ids, order, n_v_zones);
1306
1307 /* Now loop on zones in id order */
1308
1309 cs_tree_node_t *tn_head = NULL;
1310 cs_tree_node_t *tn_tail = NULL;
1311
1312 for (int i = 0; i < n_v_zones; i++) {
1313
1314 int z_id = z_ids[order[i]];
1315
1316 cs_tree_node_t *tn = _v_zone_node_by_id(tn_parent, z_id);
1317 _tree_node_remove(tn);
1318 if (tn_head == NULL) {
1319 tn_head = tn;
1320 tn_tail = tn;
1321 }
1322 else {
1323 tn->prev = tn_tail;
1324 tn_tail->next = tn;
1325 tn_tail = tn;
1326 }
1327
1328 }
1329
1330 if (tn_parent->children != NULL)
1331 tn_parent->children->prev = tn_tail;
1332 tn_tail->next = tn_parent->children;
1333 tn_parent->children = tn_head;
1334
1335 BFT_FREE(order);
1336 BFT_FREE(z_ids);
1337 }
1338
1339 /* Boundary zones */
1340 /*--------------- */
1341
1342 /* Loop on boundary condition zones */
1343
1344 tn_parent = cs_tree_get_node(cs_glob_tree,
1345 "boundary_conditions");
1346
1347 need_reorder = false;
1348 int z_id_max = 0;
1349
1350 id = 0;
1351 for (cs_tree_node_t *tn = cs_tree_node_get_child(tn_parent, "boundary");
1352 tn != NULL;
1353 tn = cs_tree_node_get_next_of_name(tn), id++) {
1354
1355 /* Zone id in tree; note that the name tag for boundary zones actually
1356 defines an integer (1 to n). This tag should be removed in the
1357 future to only use the zone label (the actual zone name). */
1358
1359 const char *id_s = cs_tree_node_get_tag(tn, "name");
1360 if (id_s != NULL) {
1361 int z_t_id = atoi(id_s);
1362 if (z_t_id != id + 1)
1363 need_reorder = true;
1364 z_id_max = CS_MAX(z_id_max, z_t_id);
1365 }
1366
1367 }
1368
1369 const int n_b_zones = id;
1370
1371 if (need_reorder) {
1372
1373 cs_lnum_t *order = NULL, *z_ids = NULL;
1374 cs_tree_node_t **tn_bcs = NULL;
1375
1376 /* Build ordering array */
1377
1378 BFT_MALLOC(z_ids, n_b_zones, cs_lnum_t);
1379 BFT_MALLOC(order, n_b_zones, cs_lnum_t);
1380 BFT_MALLOC(tn_bcs, n_b_zones, cs_tree_node_t *);
1381
1382 /* Loop on volume condition zones */
1383
1384 id = 0;
1385 for (cs_tree_node_t *tn = cs_tree_node_get_child(tn_parent, "boundary");
1386 tn != NULL;
1387 tn = cs_tree_node_get_next_of_name(tn), id++) {
1388 const char *id_s = cs_tree_node_get_tag(tn, "name");
1389 if (id_s != NULL) {
1390 z_ids[id] = atoi(id_s);
1391 }
1392 else
1393 z_ids[id] = z_id_max + 1 + id;
1394 tn_bcs[id] = tn;
1395 }
1396
1397 cs_order_lnum_allocated(NULL, z_ids, order, n_b_zones);
1398
1399 BFT_FREE(z_ids);
1400
1401 /* Now loop on zones in id order */
1402
1403 cs_tree_node_t *tn_head = NULL;
1404 cs_tree_node_t *tn_tail = NULL;
1405
1406 for (int i = 0; i < n_b_zones; i++) {
1407
1408 cs_tree_node_t *tn = tn_bcs[order[i]];
1409 _tree_node_remove(tn);
1410 if (tn_head == NULL) {
1411 tn_head = tn;
1412 tn_tail = tn;
1413 }
1414 else {
1415 tn->prev = tn_tail;
1416 tn_tail->next = tn;
1417 tn_tail = tn;
1418 }
1419
1420 }
1421
1422 if (tn_parent->children != NULL)
1423 tn_parent->children->prev = tn_tail;
1424 tn_tail->next = tn_parent->children;
1425 tn_parent->children = tn_head;
1426
1427 BFT_FREE(order);
1428 BFT_FREE(tn_bcs);
1429 }
1430 }
1431
1432 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1433
1434 /*============================================================================
1435 * Public Fortran function definitions
1436 *============================================================================*/
1437
1438 /*----------------------------------------------------------------------------
1439 * Thermal model.
1440 *
1441 * Fortran Interface:
1442 *
1443 * subroutine csther ()
1444 * *****************
1445 *
1446 *----------------------------------------------------------------------------*/
1447
CS_PROCF(csther,CSTHER)1448 void CS_PROCF (csther, CSTHER) (void)
1449 {
1450 cs_thermal_model_t *thermal = cs_get_glob_thermal_model();
1451
1452 switch(cs_gui_thermal_model()) {
1453 case 10:
1454 thermal->itherm = CS_THERMAL_MODEL_TEMPERATURE;
1455 thermal->itpscl = CS_TEMPERATURE_SCALE_CELSIUS;
1456 break;
1457 case 11:
1458 thermal->itherm = CS_THERMAL_MODEL_TEMPERATURE;
1459 thermal->itpscl = CS_TEMPERATURE_SCALE_KELVIN;
1460 break;
1461 case 12:
1462 thermal->itherm = CS_THERMAL_MODEL_TEMPERATURE;
1463 thermal->itpscl = CS_TEMPERATURE_SCALE_CELSIUS;
1464 break;
1465 case 13:
1466 thermal->itherm = CS_THERMAL_MODEL_TEMPERATURE;
1467 thermal->itpscl = CS_TEMPERATURE_SCALE_CELSIUS;
1468 break;
1469 case 20:
1470 thermal->itherm = CS_THERMAL_MODEL_ENTHALPY;
1471 thermal->itpscl = CS_TEMPERATURE_SCALE_KELVIN;
1472 break;
1473 case 30:
1474 thermal->itherm = CS_THERMAL_MODEL_TOTAL_ENERGY;
1475 thermal->itpscl = CS_TEMPERATURE_SCALE_KELVIN;
1476 break;
1477 default:
1478 thermal->itherm = CS_THERMAL_MODEL_NONE;
1479 thermal->itpscl = CS_TEMPERATURE_SCALE_NONE;
1480 break;
1481 }
1482 }
1483
1484 /*----------------------------------------------------------------------------
1485 * Turbulence model
1486 *----------------------------------------------------------------------------*/
1487
cs_gui_turb_model(void)1488 void cs_gui_turb_model(void)
1489 {
1490 cs_tree_node_t *tn_t = cs_tree_get_node(cs_glob_tree,
1491 "thermophysical_models/turbulence");
1492
1493 const char *model = cs_tree_node_get_tag(tn_t, "model");
1494 if (model == NULL)
1495 return;
1496
1497 int iwallf = -1;
1498 cs_turb_model_t *turb_mdl = cs_get_glob_turb_model();
1499 cs_turb_rans_model_t *rans_mdl = cs_get_glob_turb_rans_model();
1500
1501 if (cs_gui_strcmp(model, "off"))
1502 turb_mdl->iturb = CS_TURB_NONE;
1503 else if (cs_gui_strcmp(model, "mixing_length")) {
1504 turb_mdl->iturb = CS_TURB_MIXING_LENGTH;
1505 cs_gui_node_get_child_real(tn_t, "mixing_length_scale", &(rans_mdl->xlomlg));
1506 }
1507 else if (cs_gui_strcmp(model, "k-epsilon")) {
1508 turb_mdl->iturb = CS_TURB_K_EPSILON;
1509 cs_gui_node_get_child_int(tn_t, "wall_function", &iwallf);
1510 cs_gui_node_get_child_status_int(tn_t, "gravity_terms", &(rans_mdl->igrake));
1511 }
1512 else if (cs_gui_strcmp(model, "k-epsilon-PL")) {
1513 turb_mdl->iturb = CS_TURB_K_EPSILON_LIN_PROD;
1514 cs_gui_node_get_child_int(tn_t, "wall_function", &iwallf);
1515 cs_gui_node_get_child_status_int(tn_t, "gravity_terms", &(rans_mdl->igrake));
1516 }
1517 else if (cs_gui_strcmp(model, "Rij-epsilon")) {
1518 turb_mdl->iturb = CS_TURB_RIJ_EPSILON_LRR;
1519 cs_gui_node_get_child_int(tn_t, "wall_function", &iwallf);
1520 cs_gui_node_get_child_status_int(tn_t, "gravity_terms", &(rans_mdl->igrari));
1521 }
1522 else if (cs_gui_strcmp(model, "Rij-SSG")) {
1523 turb_mdl->iturb = CS_TURB_RIJ_EPSILON_SSG;
1524 cs_gui_node_get_child_int(tn_t, "wall_function", &iwallf);
1525 cs_gui_node_get_child_status_int(tn_t, "gravity_terms", &(rans_mdl->igrari));
1526 }
1527 else if (cs_gui_strcmp(model, "Rij-EBRSM")) {
1528 turb_mdl->iturb = CS_TURB_RIJ_EPSILON_EBRSM;
1529 cs_gui_node_get_child_int(tn_t, "wall_function", &iwallf);
1530 cs_gui_node_get_child_status_int(tn_t, "gravity_terms", &(rans_mdl->igrari));
1531 }
1532 else if (cs_gui_strcmp(model, "LES_Smagorinsky")) {
1533 turb_mdl->iturb = CS_TURB_LES_SMAGO_CONST;
1534 }
1535 else if (cs_gui_strcmp(model, "LES_dynamique")) {
1536 turb_mdl->iturb = CS_TURB_LES_SMAGO_DYN;
1537 }
1538 else if (cs_gui_strcmp(model, "LES_WALE")) {
1539 turb_mdl->iturb = CS_TURB_LES_WALE;
1540 }
1541 else if (cs_gui_strcmp(model, "v2f-phi")) {
1542 turb_mdl->iturb = CS_TURB_V2F_PHI;
1543 cs_gui_node_get_child_int(tn_t, "wall_function", &iwallf);
1544 cs_gui_node_get_child_status_int(tn_t, "gravity_terms", &(rans_mdl->igrake));
1545 }
1546 else if (cs_gui_strcmp(model, "v2f-BL-v2/k")) {
1547 turb_mdl->iturb = CS_TURB_V2F_BL_V2K;
1548 cs_gui_node_get_child_int(tn_t, "wall_function", &iwallf);
1549 cs_gui_node_get_child_status_int(tn_t, "gravity_terms", &(rans_mdl->igrake));
1550 }
1551 else if (cs_gui_strcmp(model, "k-omega-SST")) {
1552 turb_mdl->iturb = CS_TURB_K_OMEGA;
1553 cs_gui_node_get_child_int(tn_t, "wall_function", &iwallf);
1554 cs_gui_node_get_child_status_int(tn_t, "gravity_terms", &(rans_mdl->igrake));
1555 }
1556 else if (cs_gui_strcmp(model, "Spalart-Allmaras")) {
1557 turb_mdl->iturb = CS_TURB_SPALART_ALLMARAS;
1558 }
1559 else
1560 bft_error(__FILE__, __LINE__, 0,
1561 _("Invalid turbulence model: %s.\n"), model);
1562
1563 if (iwallf != -1) {
1564 cs_wall_functions_t *wall_fnt = cs_get_glob_wall_functions();
1565 wall_fnt->iwallf = (cs_wall_f_type_t)iwallf;
1566 }
1567
1568 if ( turb_mdl->iturb >= CS_TURB_RIJ_EPSILON_LRR
1569 && turb_mdl->iturb <= CS_TURB_RIJ_EPSILON_EBRSM) {
1570 const char *s
1571 = cs_tree_node_get_child_value_str(tn_t, "turbulent_diffusion_model");
1572 if (s != NULL) {
1573 if (cs_gui_strcmp(s, "shir"))
1574 rans_mdl->idirsm = 0;
1575 else if (cs_gui_strcmp(s, "daly_harlow"))
1576 rans_mdl->idirsm = 1;
1577 }
1578 }
1579
1580 #if _XML_DEBUG_
1581 bft_printf("==> %s\n", __func__);
1582 bft_printf("--model: %s\n", model);
1583 bft_printf("--iturb = %i\n", turb_mdl->iturb);
1584 bft_printf("--igrake = %i\n", rans_mdl->igrake);
1585 bft_printf("--igrari = %i\n", rans_mdl->igrari);
1586 bft_printf("--iwallf = %i\n", wall_fnt->iwallf);
1587 bft_printf("--xlomlg = %f\n", rans_mdl->xlomlg);
1588 bft_printf("--idirsm = %f\n", rans_mdl->idirsm);
1589 #endif
1590 }
1591
1592 /*----------------------------------------------------------------------------
1593 * Define reference length and reference velocity for initialization of the
1594 * turbulence variables
1595 *----------------------------------------------------------------------------*/
1596
cs_gui_turb_ref_values(void)1597 void cs_gui_turb_ref_values(void)
1598 {
1599 cs_tree_node_t *tn_t = cs_tree_get_node(cs_glob_tree,
1600 "thermophysical_models/turbulence");
1601
1602 cs_turb_model_t *turb_mdl = cs_get_glob_turb_model();
1603
1604 if (turb_mdl->iturb != 0) {
1605 const char* length_choice = NULL;
1606 cs_turb_ref_values_t *ref_values = cs_get_glob_turb_ref_values();
1607
1608 ref_values->uref = 1.; /* default if not specified */
1609
1610 cs_gui_node_get_child_real(tn_t,
1611 "reference_velocity",
1612 &(ref_values->uref));
1613
1614 length_choice = _reference_length_initialization_choice();
1615
1616 if (length_choice != NULL) {
1617 if (cs_gui_strcmp(length_choice, "prescribed"))
1618 cs_gui_node_get_child_real(tn_t,
1619 "reference_length",
1620 &(ref_values->almax));
1621 }
1622 }
1623
1624 #if _XML_DEBUG_
1625 bft_printf("==> %s\n", __func__);
1626 bft_printf("--almax = %f\n", ref_values->almax);
1627 bft_printf("--uref = %f\n", ref_values->uref);
1628 #endif
1629 }
1630
1631 /*----------------------------------------------------------------------------
1632 * Specific heat variable or constant indicator.
1633 *
1634 * Fortran Interface:
1635 *
1636 * SUBROUTINE CSCPVA (ICP)
1637 * *****************
1638 *
1639 * INTEGER ICP --> specific heat variable or constant indicator
1640 *----------------------------------------------------------------------------*/
1641
CS_PROCF(cscpva,CSCPVA)1642 void CS_PROCF (cscpva, CSCPVA) (void)
1643 {
1644 int choice;
1645 cs_fluid_properties_t *phys_pp = cs_get_glob_fluid_properties();
1646
1647 if (_properties_choice_id("specific_heat", &choice))
1648 phys_pp->icp = (choice > 0) ? 0 : -1;
1649
1650 #if _XML_DEBUG_
1651 bft_printf("==> %s\n", __func__);
1652 bft_printf("--icp = %i\n", phys_pp->icp);
1653 #endif
1654 }
1655
1656 /*----------------------------------------------------------------------------
1657 * Volumic viscosity variable or constant indicator.
1658 *
1659 * Fortran Interface:
1660 *
1661 * SUBROUTINE CSCVVVA (IVISCV)
1662 * *****************
1663 *
1664 * INTEGER IVISCV --> volumic viscosity variable or constant indicator
1665 *----------------------------------------------------------------------------*/
1666
CS_PROCF(csvvva,CSVVVA)1667 void CS_PROCF (csvvva, CSVVVA) (int *iviscv)
1668 {
1669 int choice;
1670
1671 if (_properties_choice_id("volume_viscosity", &choice))
1672 *iviscv = (choice > 0) ? 0 : -1;
1673
1674 #if _XML_DEBUG_
1675 bft_printf("==> %s\n", __func__);
1676 bft_printf("--iviscv = %i\n", *iviscv);
1677 #endif
1678 }
1679
1680 /*----------------------------------------------------------------------------
1681 * Constant or variable indicator for the user scalar molecular diffusivity
1682 *
1683 * Fortran Interface:
1684 *
1685 * subroutine csivis
1686 * *****************
1687 *----------------------------------------------------------------------------*/
1688
CS_PROCF(csivis,CSIVIS)1689 void CS_PROCF (csivis, CSIVIS) (void)
1690 {
1691 int choice1, choice2;
1692 int test1, test2;
1693
1694 const int keysca = cs_field_key_id("scalar_id");
1695 const int kivisl = cs_field_key_id("diffusivity_id");
1696 const int kscavr = cs_field_key_id("first_moment_id");
1697 const int n_fields = cs_field_n_fields();
1698
1699 cs_field_t *tf = cs_thermal_model_field();
1700
1701 if ( cs_glob_physical_model_flag[CS_PHYSICAL_MODEL_FLAG] <= 0
1702 && tf != NULL) {
1703 test1 = _properties_choice_id("thermal_conductivity", &choice1);
1704 test2 = _properties_choice_id("specific_heat", &choice2);
1705
1706 if (test1 && test2) {
1707
1708 for (int f_id = 0; f_id < n_fields; f_id++) {
1709 cs_field_t *f = cs_field_by_id(f_id);
1710 if (f->type & CS_FIELD_VARIABLE) {
1711 if (f == tf) {
1712 if (choice1 || choice2)
1713 cs_field_set_key_int(f, kivisl, 0);
1714 else
1715 cs_field_set_key_int(f, kivisl, -1);
1716 }
1717 }
1718 }
1719 }
1720 }
1721
1722 for (int f_id = 0; f_id < n_fields; f_id++) {
1723 cs_field_t *f = cs_field_by_id(f_id);
1724
1725 if ( (f->type & CS_FIELD_VARIABLE)
1726 && (f->type & CS_FIELD_USER)
1727 && (f != tf)) {
1728 int iscal = cs_field_get_key_int(f, keysca);
1729 if (iscal > 0) {
1730 if (cs_field_get_key_int(f, kscavr) < 0) {
1731 if (_scalar_properties_choice(f->name, &choice1))
1732 cs_field_set_key_int(f, kivisl, choice1 - 1);
1733 // for groundwater we impose variable property
1734 if (cs_glob_physical_model_flag[CS_GROUNDWATER] > -1)
1735 cs_field_set_key_int(f, kivisl, 0);
1736 }
1737 }
1738 }
1739 }
1740
1741 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
1742 int d_f_id = -1;
1743 /*FIXME:CK*/
1744 const char *prop_choice = _properties_choice("thermal_conductivity", NULL);
1745 if (cs_gui_strcmp(prop_choice, "user_law") ||
1746 cs_gui_strcmp(prop_choice, "predefined_law"))
1747 d_f_id = 0;
1748 cs_field_t *c_temp = cs_field_by_name("temperature");
1749 cs_field_set_key_int(c_temp, kivisl, d_f_id);
1750 }
1751 }
1752
1753 /*----------------------------------------------------------------------------
1754 * Time passing parameter.
1755 *
1756 * Fortran Interface:
1757 *
1758 * SUBROUTINE CSIDTV ()
1759 * *****************
1760 *----------------------------------------------------------------------------*/
1761
CS_PROCF(csidtv,CSIDTV)1762 void CS_PROCF (csidtv, CSIDTV) (void)
1763 {
1764 cs_time_step_options_t *time_opt = cs_get_glob_time_step_options();
1765
1766 cs_tree_node_t *tn
1767 = cs_tree_get_node(cs_glob_tree, "analysis_control/time_parameters");
1768 cs_gui_node_get_child_int(tn, "time_passing", &time_opt->idtvar);
1769
1770 #if _XML_DEBUG_
1771 bft_printf("==> %s\n", __func__);
1772 bft_printf("--idtvar = %i\n", time_opt->idtvar);
1773 #endif
1774 }
1775
1776 /*----------------------------------------------------------------------------
1777 * Hydrostatic pressure parameter.
1778 *
1779 * Fortran Interface:
1780 *
1781 * SUBROUTINE CSIPHY ()
1782 * *****************
1783 *
1784 *----------------------------------------------------------------------------*/
1785
CS_PROCF(csiphy,CSIPHY)1786 void CS_PROCF (csiphy, CSIPHY) (void)
1787 {
1788 cs_velocity_pressure_param_t *vp_param
1789 = cs_get_glob_velocity_pressure_param();
1790 int result = vp_param->iphydr;
1791 cs_tree_node_t *tn
1792 = cs_tree_find_node(cs_glob_tree,
1793 "numerical_parameters/hydrostatic_pressure");
1794 cs_gui_node_get_status_int(tn, &result);
1795 vp_param->iphydr = result;
1796
1797 #if _XML_DEBUG_
1798 bft_printf("==> %s\n", __func__);
1799 bft_printf("--iphydr = %i\n", vp_model->iphydr);
1800 #endif
1801 }
1802
1803 /*----------------------------------------------------------------------------
1804 * Hydrostatic equilibrium parameter.
1805 *
1806 * Fortran Interface:
1807 *
1808 * SUBROUTINE CSCFGP (icfgrp)
1809 * *****************
1810 *
1811 * INTEGER icfgrp --> hydrostatic equilibrium
1812 *----------------------------------------------------------------------------*/
1813
CS_PROCF(cscfgp,CSCFGP)1814 void CS_PROCF (cscfgp, CSCFGP) (int *icfgrp)
1815 {
1816 int result = *icfgrp;
1817 cs_tree_node_t *tn
1818 = cs_tree_find_node(cs_glob_tree,
1819 "numerical_parameters/hydrostatic_equilibrium/");
1820 cs_gui_node_get_status_int(tn, &result);
1821 *icfgrp = result;
1822
1823 #if _XML_DEBUG_
1824 bft_printf("==> %s\n", __func__);
1825 bft_printf("--icfgrp = %i\n", *icfgrp);
1826 #endif
1827 }
1828
1829 /*----------------------------------------------------------------------------
1830 * Restart parameters.
1831 *
1832 * Fortran Interface:
1833 *
1834 * SUBROUTINE CSISUI (NTSUIT, ILEAUX, ICCVFG)
1835 * *****************
1836 *
1837 * INTEGER NTSUIT --> checkpoint frequency
1838 * INTEGER ICCFVG --> restart with frozen field
1839 *----------------------------------------------------------------------------*/
1840
CS_PROCF(csisui,CSISUI)1841 void CS_PROCF (csisui, CSISUI) (int *ntsuit,
1842 int *iccvfg)
1843 {
1844 cs_tree_node_t *tn = cs_tree_get_node(cs_glob_tree,
1845 "calculation_management/start_restart");
1846
1847 cs_gui_node_get_child_int(tn, "restart_rescue", ntsuit);
1848
1849 cs_gui_node_get_child_status_int
1850 (tn, "restart_with_auxiliary",
1851 &(cs_glob_restart_auxiliary->read_auxiliary));
1852
1853 cs_gui_node_get_child_status_int(tn, "frozen_field", iccvfg);
1854
1855 #if _XML_DEBUG_
1856 bft_printf("==> %s\n", __func__);
1857 bft_printf("--ntsuit = %i\n", *ntsuit);
1858 bft_printf("--ileaux = %i\n", cs_glob_restart_auxiliary->read_auxiliary);
1859 bft_printf("--iccvfg = %i\n", *iccvfg);
1860 #endif
1861 }
1862
1863 /*----------------------------------------------------------------------------
1864 * Time passing parameters.
1865 *
1866 * Fortran Interface:
1867 *
1868 * SUBROUTINE CSTIME ()
1869 *
1870 *----------------------------------------------------------------------------*/
1871
CS_PROCF(cstime,CSTIME)1872 void CS_PROCF (cstime, CSTIME) (void)
1873 {
1874 /* Default, forbidden values for time step factor */
1875 double cdtmin = -1., cdtmax = -1.;
1876
1877 cs_tree_node_t *tn
1878 = cs_tree_get_node(cs_glob_tree, "analysis_control/time_parameters");
1879
1880 cs_time_step_options_t *time_opt = cs_get_glob_time_step_options();
1881 cs_time_step_t *time_stp = cs_get_glob_time_step();
1882
1883 cs_gui_node_get_child_real(tn, "time_step_ref", &(time_stp->dt_ref));
1884 cs_gui_node_get_child_real(tn, "time_step_min_factor", &cdtmin);
1885 cs_gui_node_get_child_real(tn, "time_step_max_factor", &cdtmax);
1886 cs_gui_node_get_child_real(tn, "max_courant_num", &(time_opt->coumax));
1887 cs_gui_node_get_child_real(tn, "max_fourier_num", &(time_opt->foumax));
1888 cs_gui_node_get_child_real(tn, "time_step_var", &(time_opt->varrdt));
1889 cs_gui_node_get_child_real(tn, "relaxation_coefficient", &(time_opt->relxst));
1890
1891 if (cdtmin > 0)
1892 time_opt->dtmin = cdtmin * time_stp->dt_ref;
1893 if (cdtmax > 0)
1894 time_opt->dtmax = cdtmax * time_stp->dt_ref;
1895
1896 /* We keep these two lines in case we read an old XML file... */
1897 cs_gui_node_get_child_real(tn, "time_step_min", &(time_opt->dtmin));
1898 cs_gui_node_get_child_real(tn, "time_step_max", &(time_opt->dtmax));
1899
1900 /* Stop criterion */
1901
1902 cs_real_t _t_max = -1;
1903
1904 cs_gui_node_get_child_real(tn, "maximum_time", &_t_max);
1905 if (_t_max >= 0)
1906 time_stp->t_max = _t_max;
1907 else {
1908 cs_gui_node_get_child_real(tn, "maximum_time_add", &_t_max);
1909 if (_t_max >= 0)
1910 time_stp->t_max = time_stp->t_prev + _t_max;
1911 }
1912
1913 if (_t_max < 0) {
1914 int _nt_max = -1;
1915 cs_gui_node_get_child_int(tn, "iterations", &_nt_max);
1916 if (_nt_max > -1)
1917 time_stp->nt_max = _nt_max;
1918 else {
1919 cs_gui_node_get_child_int(tn, "iterations_add", &_nt_max);
1920 if (_nt_max > -1)
1921 time_stp->nt_max = time_stp->nt_prev + _nt_max;
1922 }
1923 }
1924
1925 cs_gui_node_get_child_status_int(tn,
1926 "thermal_time_step",
1927 &(time_opt->iptlro));
1928
1929 #if _XML_DEBUG_
1930 bft_printf("==> %s\n", __func__);
1931 bft_printf("--idtvar = %i\n", time_opt->idtvar);
1932 bft_printf("--iptlro = %i\n", time_opt->iptlro);
1933 bft_printf("--ntmabs = %i\n", time_stp->nt_max);
1934 bft_printf("--dtref = %f\n", time_opt->dtref);
1935 bft_printf("--dtmin = %f\n", time_opt->dtmin);
1936 bft_printf("--dtmax = %f\n", time_opt->dtmax);
1937 bft_printf("--coumax = %f\n", time_opt->coumax);
1938 bft_printf("--foumax = %f\n", time_opt->foumax);
1939 bft_printf("--varrdt = %f\n", time_opt->varrdt);
1940 bft_printf("--relxst = %f\n", time_opt->relxst);
1941 #endif
1942 }
1943
1944 /*----------------------------------------------------------------------------
1945 * Space scheme options, linear solver precision and time step factor
1946 *----------------------------------------------------------------------------*/
1947
CS_PROCF(uinum1,UINUM1)1948 void CS_PROCF (uinum1, UINUM1) (double *cdtvar)
1949 {
1950 const int key_cal_opt_id = cs_field_key_id("var_cal_opt");
1951 const int var_key_id = cs_field_key_id("variable_id");
1952 const int keysca = cs_field_key_id("scalar_id");
1953 cs_var_cal_opt_t var_cal_opt;
1954
1955 int n_fields = cs_field_n_fields();
1956 for (int f_id = 0; f_id < n_fields; f_id++) {
1957 cs_field_t *f = cs_field_by_id(f_id);
1958 if (f->type & CS_FIELD_VARIABLE) {
1959 int j = cs_field_get_key_int(f, var_key_id) -1;
1960 cs_field_get_key_struct(f, key_cal_opt_id, &var_cal_opt);
1961
1962 const char *ref_name = f->name;
1963 if ( cs_gui_strcmp(f->name, "r11")
1964 || cs_gui_strcmp(f->name, "r22")
1965 || cs_gui_strcmp(f->name, "r33")
1966 || cs_gui_strcmp(f->name, "r12")
1967 || cs_gui_strcmp(f->name, "r23")
1968 || cs_gui_strcmp(f->name, "r13"))
1969 ref_name = "rij";
1970
1971 cs_tree_node_t *tn_v = _find_node_variable(ref_name);
1972
1973 cs_gui_node_get_child_real(tn_v, "solver_precision", &var_cal_opt.epsilo);
1974 cs_gui_node_get_child_status_int(tn_v, "flux_reconstruction",
1975 &var_cal_opt.ircflu);
1976 cs_gui_node_get_child_int(tn_v, "rhs_reconstruction",
1977 &var_cal_opt.nswrsm);
1978 cs_gui_node_get_child_int(tn_v, "verbosity", &var_cal_opt.verbosity);
1979
1980 /* For CDO equation, if non-automatic value ie != -1 */
1981 cs_equation_param_t *eqp = cs_equation_param_by_name(f->name);
1982 if ( eqp != NULL && cs_gui_is_equal_real(var_cal_opt.epsilo, -1) == 0
1983 && eqp->sles_param != NULL)
1984 eqp->sles_param->eps = var_cal_opt.epsilo;
1985
1986 /* convection scheme options */
1987 if (var_cal_opt.iconv > 0) {
1988 cs_gui_node_get_child_real(tn_v, "blending_factor",
1989 &var_cal_opt.blencv);
1990 _order_scheme_value(tn_v, &var_cal_opt.ischcv);
1991 _slope_test_value(tn_v, &var_cal_opt.isstpc);
1992 }
1993
1994 /* set field calculation options */
1995 cs_field_set_key_struct(f, key_cal_opt_id, &var_cal_opt);
1996
1997 /* only for additional variables (user or model) */
1998 int isca = cs_field_get_key_int(f, keysca);
1999 if (isca > 0) {
2000 /* time step factor */
2001 cs_gui_node_get_child_real(tn_v, "time_step_factor", &cdtvar[j]);
2002 }
2003 }
2004 }
2005
2006 #if _XML_DEBUG_
2007 bft_printf("==> %s\n", __func__);
2008 for (int f_id = 0; f_id < n_fields; f_id++) {
2009 const cs_field_t *f = cs_field_by_id(f_id);
2010 if (f->type & CS_FIELD_VARIABLE) {
2011 j = cs_field_get_key_int(f, var_key_id) -1;
2012 bft_printf("-->variable[%i] = %s\n", j, f->name);
2013 bft_printf("--blencv = %f\n", var_cal_opt.blencv);
2014 bft_printf("--epsilo = %g\n", var_cal_opt.epsilo);
2015 bft_printf("--cdtvar = %g\n", cdtvar[j]);
2016 bft_printf("--ischcv = %i\n", var_cal_opt.ischcv);
2017 bft_printf("--isstpc = %i\n", var_cal_opt.isstpc);
2018 bft_printf("--ircflu = %i\n", var_cal_opt.ircflu);
2019 bft_printf("--nswrsm = %i\n", var_cal_opt.nswrsm);
2020 }
2021 }
2022 #endif
2023 }
2024
2025 /*----------------------------------------------------------------------------
2026 * Global numerical parameters.
2027 *
2028 * Fortran Interface:
2029 *
2030 * SUBROUTINE CSNUM2 (RELAXP, IMRGRA)
2031 * *****************
2032 * DOUBLE PRECISION RELAXP --> pressure relaxation
2033 * INTEGER IMRGRA --> gradient reconstruction
2034 *----------------------------------------------------------------------------*/
2035
CS_PROCF(csnum2,CSNUM2)2036 void CS_PROCF (csnum2, CSNUM2)(double *relaxp,
2037 int *imrgra)
2038 {
2039 cs_velocity_pressure_param_t *vp_param = cs_get_glob_velocity_pressure_param();
2040 cs_velocity_pressure_model_t *vp_model = cs_get_glob_velocity_pressure_model();
2041
2042 const char *choice = NULL;
2043
2044 cs_tree_node_t *tn_n = cs_tree_get_node(cs_glob_tree, "numerical_parameters");
2045
2046 int _imrgra = -1;
2047
2048 cs_ext_neighborhood_type_t enh_type = cs_ext_neighborhood_get_type();
2049
2050 choice = cs_tree_node_get_tag(cs_tree_get_node(tn_n,
2051 "gradient_reconstruction"),
2052 "choice");
2053 if (cs_gui_strcmp(choice, "green_iter"))
2054 _imrgra = 0;
2055 else if (cs_gui_strcmp(choice, "lsq"))
2056 _imrgra = 1;
2057 else if (cs_gui_strcmp(choice, "green_lsq"))
2058 _imrgra = 4;
2059 else if (cs_gui_strcmp(choice, "green_vtx"))
2060 _imrgra = 7;
2061
2062 if (_imrgra != 0 && _imrgra != 7) {
2063 choice = cs_tree_node_get_tag(cs_tree_get_node(tn_n,
2064 "extended_neighborhood"),
2065 "choice");
2066 if (cs_gui_strcmp(choice, "none")) {
2067 enh_type = CS_EXT_NEIGHBORHOOD_NONE;
2068 }
2069 else if (cs_gui_strcmp(choice, "complete")) {
2070 enh_type = CS_EXT_NEIGHBORHOOD_COMPLETE;
2071 _imrgra += 1;
2072 }
2073 else if (cs_gui_strcmp(choice, "cell_center_opposite")) {
2074 enh_type = CS_EXT_NEIGHBORHOOD_CELL_CENTER_OPPOSITE;
2075 _imrgra += 2;
2076 }
2077 else if (cs_gui_strcmp(choice, "non_ortho_max")) {
2078 enh_type = CS_EXT_NEIGHBORHOOD_NON_ORTHO_MAX;
2079 _imrgra += 2;
2080 }
2081 }
2082
2083 cs_ext_neighborhood_set_type(enh_type);
2084
2085 if (_imrgra > -1)
2086 *imrgra = _imrgra;
2087
2088 int _idilat = -1;
2089
2090 choice = cs_tree_node_get_tag(cs_tree_get_node(tn_n, "algo_density_variation"),
2091 "choice");
2092 if (cs_gui_strcmp(choice, "boussi"))
2093 _idilat = 0;
2094 else if (cs_gui_strcmp(choice, "dilat_std"))
2095 _idilat = 1;
2096 else if (cs_gui_strcmp(choice, "dilat_unstd"))
2097 _idilat = 2;
2098 else if (cs_gui_strcmp(choice, "low_mach"))
2099 _idilat = 3;
2100 else if (cs_gui_strcmp(choice, "algo_fire"))
2101 _idilat = 4;
2102
2103 if (_idilat > -1)
2104 vp_model->idilat = _idilat;
2105
2106 _numerical_int_parameters("gradient_transposed", &(vp_model->ivisse));
2107 _numerical_int_parameters("velocity_pressure_coupling", &(vp_param->ipucou));
2108 _numerical_int_parameters("piso_sweep_number", &(vp_param->nterup));
2109 _numerical_double_parameters("pressure_relaxation", relaxp);
2110
2111 #if _XML_DEBUG_
2112 bft_printf("==> %s\n", __func__);
2113 bft_printf("--ivisse = %i\n", vp_model->ivisse);
2114 bft_printf("--ipucou = %i\n", vp_model->ipucou);
2115 bft_printf("--imrgra = %i\n", *imrgra);
2116 bft_printf("--nterup = %i\n", vp_param->nterup);
2117 bft_printf("--relaxp = %f\n", *relaxp);
2118 #endif
2119 }
2120
2121 /*----------------------------------------------------------------------------
2122 * Treatment of gravity and fluid physical properties
2123 * Initialize reference pressure and temperature if present
2124 *----------------------------------------------------------------------------*/
2125
2126 void
cs_gui_physical_properties(void)2127 cs_gui_physical_properties(void)
2128 {
2129 int choice;
2130 const char *material = NULL;
2131
2132 const int itherm = cs_glob_thermal_model->itherm;
2133
2134 cs_physical_constants_t *phys_cst = cs_get_glob_physical_constants();
2135
2136 _gravity_value("gravity_x", &(phys_cst->gravity[0]));
2137 _gravity_value("gravity_y", &(phys_cst->gravity[1]));
2138 _gravity_value("gravity_z", &(phys_cst->gravity[2]));
2139
2140 cs_real_t w_x, w_y, w_z;
2141 w_x = 0.;
2142 w_y = 0.;
2143 w_z = 0.;
2144
2145 _coriolis_value("omega_x", &w_x);
2146 _coriolis_value("omega_y", &w_y);
2147 _coriolis_value("omega_z", &w_z);
2148
2149 if (w_x*w_x + w_y*w_y + w_z*w_z > 0.) {
2150 cs_rotation_define(w_x, w_y, w_z, 0, 0, 0);
2151 phys_cst->icorio = 1;
2152 }
2153 else
2154 phys_cst->icorio = 0;
2155
2156 cs_fluid_properties_t *phys_pp = cs_get_glob_fluid_properties();
2157 cs_gui_fluid_properties_value("reference_pressure", &(phys_pp->p0));
2158
2159 /* Variable rho and viscl */
2160 if (_properties_choice_id("density", &choice))
2161 phys_pp->irovar = choice;
2162
2163 if (_properties_choice_id("molecular_viscosity", &choice))
2164 phys_pp->ivivar = choice;
2165 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1)
2166 if (_properties_choice_id("molecular_viscosity", &choice))
2167 phys_pp->ivivar = choice;
2168
2169 /* Read T0 in each case for user */
2170 cs_gui_fluid_properties_value("reference_temperature", &(phys_pp->t0));
2171
2172 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1)
2173 cs_gui_fluid_properties_value("reference_molar_mass", &(phys_pp->xmasmr));
2174
2175 material = _thermal_table_choice("material");
2176 if (material != NULL) {
2177 if (!(cs_gui_strcmp(material, "user_material"))) {
2178 cs_phys_prop_thermo_plane_type_t thermal_plane = CS_PHYS_PROP_PLANE_PH;
2179 if (itherm <= CS_THERMAL_MODEL_TEMPERATURE)
2180 thermal_plane = CS_PHYS_PROP_PLANE_PT;
2181 //else if (itherm == CS_THERMAL_MODEL_TOTAL_ENERGY)
2182 // // TODO compressible
2183 // thermal_plane = CS_PHYS_PROP_PLANE_PS;
2184
2185 const int itpscl = cs_glob_thermal_model->itpscl;
2186
2187 cs_thermal_table_set(material,
2188 _thermal_table_choice("method"),
2189 _thermal_table_option("reference"),
2190 thermal_plane,
2191 itpscl);
2192 }
2193 }
2194
2195 cs_vof_parameters_t *vof_param = cs_get_glob_vof_parameters();
2196
2197 if (_thermal_table_needed("density") == 0) {
2198 cs_gui_properties_value("density", &phys_pp->ro0);
2199 if (vof_param->vof_model & CS_VOF_ENABLED) {
2200 cs_gui_properties_value_by_fluid_id(0, "density", &vof_param->rho1);
2201 cs_gui_properties_value_by_fluid_id(1, "density", &vof_param->rho2);
2202 }
2203 }
2204 else {
2205 cs_phys_prop_compute(CS_PHYS_PROP_DENSITY,
2206 1,
2207 0,
2208 0,
2209 &phys_pp->p0,
2210 &phys_pp->t0,
2211 &phys_pp->ro0);
2212 }
2213
2214 const char *mv_name = "molecular_viscosity";
2215 if (_thermal_table_needed(mv_name) == 0) {
2216 cs_gui_properties_value(mv_name, &phys_pp->viscl0);
2217 if (vof_param->vof_model & CS_VOF_ENABLED) {
2218 cs_gui_properties_value_by_fluid_id(0, mv_name, &vof_param->mu1);
2219 cs_gui_properties_value_by_fluid_id(1, mv_name, &vof_param->mu2);
2220 }
2221 }
2222 else {
2223 cs_phys_prop_compute(CS_PHYS_PROP_DYNAMIC_VISCOSITY,
2224 1,
2225 0,
2226 0,
2227 &phys_pp->p0,
2228 &phys_pp->t0,
2229 &phys_pp->viscl0);
2230 }
2231
2232 if (_thermal_table_needed("specific_heat") == 0)
2233 cs_gui_properties_value("specific_heat", &phys_pp->cp0);
2234 else
2235 cs_phys_prop_compute(CS_PHYS_PROP_ISOBARIC_HEAT_CAPACITY,
2236 1,
2237 0,
2238 0,
2239 &phys_pp->p0,
2240 &phys_pp->t0,
2241 &phys_pp->cp0);
2242
2243 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
2244 cs_gui_properties_value("volume_viscosity", &phys_pp->viscv0);
2245 double visls_0 = -1;
2246 cs_gui_properties_value("thermal_conductivity", &visls_0);
2247 cs_field_set_key_double(cs_field_by_name("temperature"),
2248 cs_field_key_id("diffusivity_ref"),
2249 visls_0);
2250 }
2251
2252 int n_zones_pp
2253 = cs_volume_zone_n_type_zones(CS_VOLUME_ZONE_PHYSICAL_PROPERTIES);
2254 if (n_zones_pp > 1) {
2255 phys_pp->irovar = 1;
2256 phys_pp->ivivar = 1;
2257
2258 cs_field_t *tf = cs_thermal_model_field();
2259 if (tf != NULL) {
2260 phys_pp->icp = 1;
2261 int k = cs_field_key_id("diffusivity_id");
2262 int cond_diff_id = cs_field_get_key_int(tf, k);
2263 if (cond_diff_id < 0)
2264 cs_field_set_key_int(tf, k, 0);
2265 }
2266 }
2267
2268 #if _XML_DEBUG_
2269 bft_printf("==> %s\n", __func__);
2270 bft_printf("--gx = %f \n", phys_cst->gravity[0]);
2271 bft_printf("--gy = %f \n", phys_cst->gravity[1]);
2272 bft_printf("--gz = %f \n", phys_cst->gravity[2]);
2273 bft_printf("--icorio = %i \n", cs_glob_physical_constants->icorio);
2274 bft_printf("--rho = %g , variable %i\n",
2275 cs_glob_fluid_properties->ro0,
2276 cs_glob_fluid_properties->irovar);
2277 bft_printf("--mu = %g , variable %i \n",
2278 cs_glob_fluid_properties->viscl0,
2279 cs_glob_fluid_properties->ivivar);
2280 bft_printf("--Cp = %g \n", cs_glob_fluid_properties->cp0);
2281 bft_printf("--T0 = %f \n", cs_glob_fluid_properties->t0);
2282 bft_printf("--P0 = %f \n", cs_glob_fluid_properties->p0);
2283 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
2284 bft_printf("--viscv0 = %g \n", *viscv0);
2285 bft_printf("--xmasmr = %f \n", cs_glob_fluid_properties->xmasmr);
2286 }
2287 #endif
2288 }
2289
2290 /*----------------------------------------------------------------------------
2291 * Read minimum / maximum values (used in clipping) and turbulent flux model
2292 * for additional user or model variables.
2293 *----------------------------------------------------------------------------*/
2294
CS_PROCF(cssca2,CSSCA2)2295 void CS_PROCF (cssca2, CSSCA2) (void)
2296 {
2297 #if _XML_DEBUG_
2298 bft_printf("==> %s\n", __func__);
2299 #endif
2300
2301 const cs_turb_model_t *turb_model = cs_get_glob_turb_model();
2302 assert(turb_model != NULL);
2303
2304 const int kscmin = cs_field_key_id("min_scalar_clipping");
2305 const int kscmax = cs_field_key_id("max_scalar_clipping");
2306
2307 /* Specific physics: the min max of the model scalar are not given */
2308 const int keysca = cs_field_key_id("scalar_id");
2309 const int kturt = cs_field_key_id("turbulent_flux_model");
2310
2311 for (int f_id = 0; f_id < cs_field_n_fields(); f_id++) {
2312 cs_field_t *f = cs_field_by_id(f_id);
2313 if (f->type & CS_FIELD_VARIABLE) { /* variable ? */
2314 int i = cs_field_get_key_int(f, keysca) - 1;
2315 if (i > -1) { /* additional user or model variable ? */
2316
2317 double scal_min = cs_field_get_key_double(f, kscmin);
2318 double scal_max = cs_field_get_key_double(f, kscmax);
2319
2320 cs_tree_node_t *tn_v = _find_node_variable(f->name);
2321 if (tn_v != NULL) { /* variable is in xml ? */
2322 cs_gui_node_get_child_real(tn_v, "min_value", &scal_min);
2323 cs_gui_node_get_child_real(tn_v, "max_value", &scal_max);
2324 cs_field_set_key_double(f, kscmin, scal_min);
2325 cs_field_set_key_double(f, kscmax, scal_max);
2326
2327 #if _XML_DEBUG_
2328 bft_printf("--min_scalar_clipping[%i] = %f\n", i, scal_min);
2329 bft_printf("--max_scalar_clipping[%i] = %f\n", i, scal_max);
2330 #endif
2331
2332 if (turb_model->order == CS_TURB_SECOND_ORDER) {
2333 int turb_mdl;
2334 _variable_turbulent_flux_model(tn_v, &turb_mdl);
2335 cs_field_set_key_int(f, kturt, turb_mdl);
2336 #if _XML_DEBUG_
2337 bft_printf("--turb_model[%i] = %d\n", i, turb_mdl);
2338 #endif
2339 }
2340
2341 }
2342 }
2343 }
2344 }
2345 }
2346
2347 /*----------------------------------------------------------------------------
2348 * Read reference dynamic and user scalar viscosity
2349 *----------------------------------------------------------------------------*/
2350
CS_PROCF(cssca3,CSSCA3)2351 void CS_PROCF (cssca3, CSSCA3) (void)
2352 {
2353 double result, density;
2354
2355 const int kscavr = cs_field_key_id("first_moment_id");
2356 const int kvisls0 = cs_field_key_id("diffusivity_ref");
2357
2358 const int itherm = cs_glob_thermal_model->itherm;
2359
2360 cs_fluid_properties_t *fprops = cs_get_glob_fluid_properties();
2361
2362 if (itherm != CS_THERMAL_MODEL_NONE) {
2363
2364 if (_thermal_table_needed("thermal_conductivity") == 0)
2365 cs_gui_properties_value("thermal_conductivity", &(fprops->lambda0));
2366 else
2367 cs_phys_prop_compute(CS_PHYS_PROP_THERMAL_CONDUCTIVITY,
2368 1,
2369 0,
2370 0,
2371 &(cs_glob_fluid_properties->p0),
2372 &(cs_glob_fluid_properties->t0),
2373 &(fprops->lambda0));
2374
2375 double visls_0 = fprops->lambda0;
2376
2377 /* for the Temperature, the diffusivity factor is not divided by Cp
2378 * i.e. it remains lambda */
2379 if (itherm != CS_THERMAL_MODEL_TEMPERATURE)
2380 visls_0 /= cs_glob_fluid_properties->cp0;
2381
2382 cs_field_t *tf = cs_thermal_model_field();
2383 cs_field_set_key_double(tf, kvisls0, visls_0);
2384
2385 }
2386
2387 /* User scalar
2388 In the interface, the user gives the diffusion coefficient, whereas in
2389 the solver, one sets the diffusivity, thus one need to multiply
2390 this coefficient by the density to remain coherent */
2391
2392 if (cs_glob_physical_model_flag[CS_GROUNDWATER] < 0) {
2393 int n_fields = cs_field_n_fields();
2394 for (int f_id = 0; f_id < n_fields; f_id++) {
2395 cs_field_t *f = cs_field_by_id(f_id);
2396 if ( (f->type & CS_FIELD_VARIABLE)
2397 && (f->type & CS_FIELD_USER)) {
2398 if (cs_field_get_key_int(f, kscavr) < 0) {
2399
2400 if ( cs_glob_physical_model_flag[CS_COMBUSTION_PCLC] > -1
2401 || cs_glob_physical_model_flag[CS_COMBUSTION_COAL] > -1) {
2402 /* Air molar mass */
2403 result = 0.028966;
2404 cs_gui_fluid_properties_value("reference_molar_mass", &result);
2405 if (result <= 0)
2406 bft_error
2407 (__FILE__, __LINE__, 0,
2408 _("mass molar value is zero or not found in the xml file.\n"));
2409 density = cs_glob_fluid_properties->p0 *
2410 result / (8.31446 *(cs_glob_fluid_properties->t0));
2411 }
2412 else
2413 density = cs_glob_fluid_properties->ro0;
2414
2415 double visls_0 = cs_field_get_key_double(f, kvisls0);
2416 double coeff = visls_0 / density;
2417 _scalar_diffusion_value(f, &coeff);
2418 visls_0 = coeff * density;
2419
2420 cs_field_set_key_double(f, kvisls0, visls_0);
2421 }
2422 }
2423 }
2424 }
2425 }
2426
2427 /*----------------------------------------------------------------------------
2428 * Define porosity.
2429 *
2430 * Fortran Interface:
2431 *
2432 * SUBROUTINE UIPORO
2433 *----------------------------------------------------------------------------*/
2434
CS_PROCF(uiporo,UIPORO)2435 void CS_PROCF(uiporo, UIPORO)(void)
2436 {
2437 const cs_lnum_t n_cells_ext = cs_glob_mesh->n_cells_with_ghosts;
2438
2439 int n_zones = cs_volume_zone_n_zones();
2440
2441 /* Porosity fields */
2442 cs_field_t *fporo = CS_F_(poro);
2443 cs_field_t *ftporo = CS_F_(t_poro);
2444
2445 if (fporo != NULL)
2446 cs_array_set_value_real(n_cells_ext, 1, 1., fporo->val);
2447
2448 if (ftporo != NULL) {
2449 cs_real_6_t *porosf = (cs_real_6_t *)ftporo->val;
2450 for (cs_lnum_t iel = 0; iel < n_cells_ext; iel++) {
2451 porosf[iel][0] = 1.;
2452 porosf[iel][1] = 1.;
2453 porosf[iel][2] = 1.;
2454 porosf[iel][3] = 0.;
2455 porosf[iel][4] = 0.;
2456 porosf[iel][5] = 0.;
2457 }
2458 }
2459
2460 cs_tree_node_t *tn_p
2461 = cs_tree_get_node(cs_glob_tree,
2462 "thermophysical_models/porosities/porosity");
2463
2464 for (int z_id = 0; z_id < n_zones; z_id++) {
2465 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
2466
2467 if (z->type & CS_VOLUME_ZONE_POROSITY) {
2468
2469 cs_tree_node_t *tn_zp = _add_zone_id_test_attribute(tn_p, z->id);
2470 const char *mdl = cs_tree_node_get_child_value_str(tn_zp, "model");
2471 const char *formula = cs_tree_node_get_child_value_str(tn_zp, "formula");
2472
2473 if (formula != NULL) {
2474 if (cs_gui_strcmp(mdl, "anisotropic")) {
2475 cs_field_t *fmeg[2] = {fporo, ftporo};
2476 cs_meg_volume_function(z, fmeg);
2477
2478 } else {
2479 cs_field_t *fmeg[1] = {fporo};
2480 cs_meg_volume_function(z, fmeg);
2481 }
2482
2483 }
2484 }
2485 }
2486
2487 cs_porous_model_auto_face_porosity();
2488 }
2489
2490 /*----------------------------------------------------------------------------
2491 * User momentum source terms.
2492 *
2493 * Fortran Interface:
2494 *
2495 * subroutine uitsnv (vel, tsexp, tsimp)
2496 * *****************
2497 *
2498 * double precision vel <-- fluid velocity
2499 * double precision tsexp --> explicit source terms
2500 * double precision tsimp --> implicit source terms
2501 *----------------------------------------------------------------------------*/
2502
CS_PROCF(uitsnv,UITSNV)2503 void CS_PROCF(uitsnv, UITSNV)(const cs_real_3_t *restrict vel,
2504 cs_real_3_t *restrict tsexp,
2505 cs_real_33_t *restrict tsimp)
2506 {
2507 const cs_real_t *restrict cell_f_vol = cs_glob_mesh_quantities->cell_f_vol;
2508
2509 double Su, Sv, Sw;
2510 double dSudu, dSudv, dSudw;
2511 double dSvdu, dSvdv, dSvdw;
2512 double dSwdu, dSwdv, dSwdw;
2513
2514 #if _XML_DEBUG_
2515 bft_printf("==> %s\n", __func__);
2516 #endif
2517
2518 int n_zones = cs_volume_zone_n_zones();
2519
2520 cs_tree_node_t *tn_mf
2521 = cs_tree_get_node(cs_glob_tree,
2522 "thermophysical_models/source_terms/momentum_formula");
2523
2524 for (int z_id = 0; z_id < n_zones; z_id++) {
2525 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
2526
2527 if (! (z->type & CS_VOLUME_ZONE_SOURCE_TERM))
2528 continue;
2529
2530 if (_zone_id_is_type(z->id, "momentum_source_term")) {
2531 const cs_lnum_t n_cells = z->n_elts;
2532 const cs_lnum_t *cell_ids = z->elt_ids;
2533
2534 cs_tree_node_t *tn = _add_zone_id_test_attribute(tn_mf, z->id);
2535 const char *formula = cs_tree_node_get_value_str(tn);
2536
2537 if (formula != NULL) {
2538
2539 cs_real_t *st_vals = cs_meg_source_terms(z,
2540 "momentum",
2541 "momentum_source_term");
2542
2543 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
2544 cs_lnum_t c_id = cell_ids[e_id];
2545
2546 /* Read values from the newly created array */
2547 Su = st_vals[12*e_id];
2548 Sv = st_vals[12*e_id + 1];
2549 Sw = st_vals[12*e_id + 2];
2550
2551 dSudu = st_vals[12*e_id + 3];
2552 dSudv = st_vals[12*e_id + 4];
2553 dSudw = st_vals[12*e_id + 5];
2554
2555 dSvdu = st_vals[12*e_id + 6];
2556 dSvdv = st_vals[12*e_id + 7];
2557 dSvdw = st_vals[12*e_id + 8];
2558
2559 dSwdu = st_vals[12*e_id + 9];
2560 dSwdv = st_vals[12*e_id + 10];
2561 dSwdw = st_vals[12*e_id + 11];
2562
2563 /* Fill the explicit and implicit source terms' arrays */
2564 tsexp[c_id][0] = cell_f_vol[c_id]
2565 * ( Su
2566 - dSudu * vel[c_id][0]
2567 - dSudv * vel[c_id][1]
2568 - dSudw * vel[c_id][2] );
2569
2570 tsexp[c_id][1] = cell_f_vol[c_id]
2571 * ( Sv
2572 - dSvdu * vel[c_id][0]
2573 - dSvdv * vel[c_id][1]
2574 - dSvdw * vel[c_id][2] );
2575
2576 tsexp[c_id][2] = cell_f_vol[c_id]
2577 * ( Sw
2578 - dSwdu * vel[c_id][0]
2579 - dSwdv * vel[c_id][1]
2580 - dSwdw * vel[c_id][2] );
2581
2582 tsimp[c_id][0][0] = cell_f_vol[c_id]*dSudu;
2583 tsimp[c_id][0][1] = cell_f_vol[c_id]*dSudv;
2584 tsimp[c_id][0][2] = cell_f_vol[c_id]*dSudw;
2585 tsimp[c_id][1][0] = cell_f_vol[c_id]*dSvdu;
2586 tsimp[c_id][1][1] = cell_f_vol[c_id]*dSvdv;
2587 tsimp[c_id][1][2] = cell_f_vol[c_id]*dSvdw;
2588 tsimp[c_id][2][0] = cell_f_vol[c_id]*dSwdu;
2589 tsimp[c_id][2][1] = cell_f_vol[c_id]*dSwdv;
2590 tsimp[c_id][2][2] = cell_f_vol[c_id]*dSwdw;
2591
2592 }
2593 if (st_vals != NULL)
2594 BFT_FREE(st_vals);
2595 }
2596 }
2597 }
2598 }
2599
2600 /*----------------------------------------------------------------------------
2601 * User scalar source terms.
2602 *
2603 * Fortran Interface:
2604 *
2605 * subroutine uitssc (f_id, pvar, tsexp, tsimp)
2606 * *****************
2607 *
2608 * integer idarcy <-- groundwater module activation
2609 * integer f_id <-- field id
2610 * double precision pvar <-- scalar
2611 * double precision tsexp --> explicit source terms
2612 * double precision tsimp --> implicit source terms
2613 *----------------------------------------------------------------------------*/
2614
CS_PROCF(uitssc,UITSSC)2615 void CS_PROCF(uitssc, UITSSC)(const int *idarcy,
2616 const int *f_id,
2617 const cs_real_t *restrict pvar,
2618 cs_real_t *restrict tsexp,
2619 cs_real_t *restrict tsimp)
2620 {
2621 const cs_real_t *restrict cell_f_vol = cs_glob_mesh_quantities->cell_f_vol;
2622
2623 const char *formula = NULL;
2624
2625 cs_field_t *f = cs_field_by_id(*f_id);
2626
2627 #if _XML_DEBUG_
2628 bft_printf("==> %s\n", __func__);
2629 #endif
2630
2631 int n_zones = cs_volume_zone_n_zones();
2632
2633 for (int z_id = 0; z_id < n_zones; z_id++) {
2634 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
2635
2636 if (! (z->type & CS_VOLUME_ZONE_SOURCE_TERM))
2637 continue;
2638
2639 /* species source term */
2640 if (_zone_id_is_type(z->id, "scalar_source_term")) {
2641 const cs_lnum_t n_cells = z->n_elts;
2642 const cs_lnum_t *cell_ids = z->elt_ids;
2643
2644 cs_tree_node_t *tn
2645 = cs_tree_get_node(cs_glob_tree,
2646 "thermophysical_models/source_terms/scalar_formula");
2647 char z_id_str[32];
2648 snprintf(z_id_str, 31, "%d", z->id);
2649 while (tn != NULL){
2650 const char *name = cs_gui_node_get_tag(tn, "name");
2651 const char *zone_id = cs_gui_node_get_tag(tn, "zone_id");
2652 if (cs_gui_strcmp(name, f->name) && cs_gui_strcmp(zone_id, z_id_str))
2653 break;
2654 tn = cs_tree_node_get_next_of_name(tn);
2655 }
2656 formula = cs_tree_node_get_value_str(tn);
2657
2658 if (formula != NULL) {
2659 cs_real_t *st_vals = cs_meg_source_terms(z,
2660 f->name,
2661 "scalar_source_term");
2662
2663 cs_real_t sign = 1.0;
2664 cs_real_t non_linear = 1.0;
2665 /* for groundwater flow, the user filled in the positive radioactive
2666 decay rate (lambda) - this source term is always linear:
2667 -lambda Y^{n+1} */
2668 if (*idarcy > -1) {
2669 sign = -1.0;
2670 non_linear = 0.;
2671 }
2672
2673 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
2674 cs_lnum_t c_id = cell_ids[e_id];
2675 tsimp[c_id] = cell_f_vol[c_id] * sign * st_vals[2 * e_id + 1];
2676 tsexp[c_id] = cell_f_vol[c_id] * st_vals[2 * e_id]
2677 - non_linear * tsimp[c_id] * pvar[c_id];
2678 }
2679 if (st_vals != NULL)
2680 BFT_FREE(st_vals);
2681 }
2682 }
2683 }
2684 }
2685
2686 /*----------------------------------------------------------------------------
2687 * Thermal scalar source terms.
2688 *
2689 * Fortran Interface:
2690 *
2691 * subroutine uitsth (f_id, pvar, tsexp, tsimp)
2692 * *****************
2693 *
2694 * integer f_id <-- field id
2695 * double precision pvar <-- scalar
2696 * double precision tsexp --> explicit source terms
2697 * double precision tsimp --> implicit source terms
2698 *----------------------------------------------------------------------------*/
2699
CS_PROCF(uitsth,UITSTH)2700 void CS_PROCF(uitsth, UITSTH)(const int *f_id,
2701 const cs_real_t *restrict pvar,
2702 cs_real_t *restrict tsexp,
2703 cs_real_t *restrict tsimp)
2704 {
2705 const cs_real_t *restrict cell_f_vol = cs_glob_mesh_quantities->cell_f_vol;
2706
2707 const char *formula = NULL;
2708
2709 cs_field_t *f = cs_field_by_id(*f_id);
2710
2711 /* number of volumic zone */
2712
2713 #if _XML_DEBUG_
2714 bft_printf("==> %s\n", __func__);
2715 #endif
2716
2717 int n_zones = cs_volume_zone_n_zones();
2718
2719 for (int z_id = 0; z_id < n_zones; z_id++) {
2720 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
2721
2722 if (!(z->type & CS_VOLUME_ZONE_SOURCE_TERM))
2723 continue;
2724
2725 /* species source term */
2726 if (_zone_id_is_type(z->id, "thermal_source_term")) {
2727 const cs_lnum_t n_cells = z->n_elts;
2728 const cs_lnum_t *cell_ids = z->elt_ids;
2729
2730 cs_tree_node_t *tn
2731 = cs_tree_get_node(cs_glob_tree,
2732 "thermophysical_models/source_terms/thermal_formula");
2733 char z_id_str[32];
2734 snprintf(z_id_str, 31, "%d", z->id);
2735 while (tn != NULL) {
2736
2737 const char *name = cs_gui_node_get_tag(tn, "name");
2738 const char *zone_id = cs_gui_node_get_tag(tn, "zone_id");
2739 if (cs_gui_strcmp(name, f->name) && cs_gui_strcmp(zone_id, z_id_str))
2740 break;
2741 tn = cs_tree_node_get_next_of_name(tn);
2742 }
2743 formula = cs_tree_node_get_value_str(tn);
2744
2745 if (formula != NULL) {
2746
2747 cs_real_t *st_vals = cs_meg_source_terms(z,
2748 f->name,
2749 "thermal_source_term");
2750
2751 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
2752 cs_lnum_t c_id = cell_ids[e_id];
2753
2754 tsimp[c_id] = cell_f_vol[c_id] * st_vals[2 * e_id + 1];
2755 tsexp[c_id] = cell_f_vol[c_id] * st_vals[2 * e_id]
2756 - tsimp[c_id] * pvar[c_id];
2757 }
2758 if (st_vals != NULL)
2759 BFT_FREE(st_vals);
2760 }
2761 }
2762 }
2763 }
2764
2765 /*----------------------------------------------------------------------------
2766 * Variables and user scalars initialization.
2767 *
2768 * Fortran Interface:
2769 *
2770 * subroutine uiiniv
2771 * *****************
2772 *
2773 * integer isuite <-- restart indicator
2774 * integer idarcy <-- groundwater module activation
2775 * integer iccfth --> type of initialization (compressible model)
2776 *----------------------------------------------------------------------------*/
2777
CS_PROCF(uiiniv,UIINIV)2778 void CS_PROCF(uiiniv, UIINIV)(const int *isuite,
2779 const int *idarcy,
2780 int *iccfth)
2781 {
2782 /* Coal combustion: the initialization of the model scalar are not given */
2783
2784 int ccfth = 0;
2785
2786 #if _XML_DEBUG_
2787 bft_printf("==> %s\n", __func__);
2788 #endif
2789
2790 const int n_zones = cs_volume_zone_n_zones();
2791
2792 for (int z_id = 0; z_id < n_zones; z_id++) {
2793 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
2794
2795 if (z->type & CS_VOLUME_ZONE_INITIALIZATION) {
2796 const cs_lnum_t n_cells = z->n_elts;
2797 const cs_lnum_t *cell_ids = z->elt_ids;
2798
2799 char z_id_str[32];
2800 snprintf(z_id_str, 31, "%d", z_id);
2801
2802 if (*isuite == 0) {
2803
2804 cs_tree_node_t *tn_velocity
2805 = cs_tree_get_node(cs_glob_tree,
2806 "thermophysical_models/velocity_pressure"
2807 "/initialization/formula");
2808 tn_velocity = _add_zone_id_test_attribute(tn_velocity, z->id);
2809 const char *formula_uvw = cs_tree_node_get_value_str(tn_velocity);
2810
2811 cs_field_t *c_vel = cs_field_by_name("velocity");
2812
2813 if (formula_uvw != NULL) {
2814 cs_real_t *ini_vals = cs_meg_initialization(z, "velocity");
2815 if (ini_vals != NULL) {
2816 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
2817 cs_lnum_t c_id = cell_ids[e_id];
2818 for (int d = 0; d < 3; d++)
2819 c_vel->val[3 * c_id + d] = ini_vals[3 * e_id + d];
2820 }
2821 BFT_FREE(ini_vals);
2822 }
2823 }
2824 else {
2825 for (cs_lnum_t icel = 0; icel < n_cells; icel++) {
2826 cs_lnum_t iel = cell_ids[icel];
2827 for (cs_lnum_t j = 0; j < 3; j++)
2828 c_vel->val[3 * iel + j] = 0.0;
2829 }
2830 }
2831
2832 /* pressure initialization for groundwater model */
2833 if (*idarcy > 0) {
2834 const char *formula = NULL;
2835
2836 cs_tree_node_t *tn = _find_node_variable("hydraulic_head");
2837 tn = cs_tree_find_node(tn, "formula");
2838 tn = _add_zone_id_test_attribute(tn, z->id);
2839 formula = cs_tree_node_get_value_str(tn);
2840
2841 cs_field_t *c = cs_field_by_name_try("hydraulic_head");
2842
2843 if (formula != NULL) {
2844 cs_real_t *ini_vals = cs_meg_initialization(z, "hydraulic_head");
2845 if (ini_vals != NULL) {
2846 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
2847 cs_lnum_t c_id = cell_ids[e_id];
2848 c->val[c_id] = ini_vals[e_id];
2849 }
2850 BFT_FREE(ini_vals);
2851 }
2852 }
2853
2854 }
2855
2856 /* Turbulence variables initialization */
2857 const char *choice = _turbulence_initialization_choice(z_id_str);
2858
2859 if (cs_gui_strcmp(choice, "formula")) {
2860
2861 const char *formula_turb = NULL;
2862 cs_tree_node_t *tn_turb
2863 = cs_tree_get_node(cs_glob_tree,
2864 "thermophysical_models/"
2865 "turbulence/initialization");
2866 tn_turb = _add_zone_id_test_attribute(tn_turb, z->id);
2867 tn_turb = cs_tree_get_node(tn_turb, "formula");
2868 formula_turb = cs_tree_node_get_value_str(tn_turb);
2869
2870 if (formula_turb != NULL) {
2871 const char *model = cs_gui_get_thermophysical_model("turbulence");
2872 if (model == NULL)
2873 break;
2874 if (cs_gui_strcmp(model, "off"))
2875 break;
2876
2877 cs_real_t *ini_vals = cs_meg_initialization(z, "turbulence");
2878
2879 if (ini_vals != NULL) {
2880
2881 if ( cs_gui_strcmp(model, "k-epsilon")
2882 || cs_gui_strcmp(model, "k-epsilon-PL")) {
2883
2884 cs_field_t *c_k = cs_field_by_name("k");
2885 cs_field_t *c_eps = cs_field_by_name("epsilon");
2886
2887 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
2888 cs_lnum_t c_id = cell_ids[e_id];
2889 c_k->val[c_id] = ini_vals[2 * e_id];
2890 c_eps->val[c_id] = ini_vals[2 * e_id + 1];
2891 }
2892 }
2893 else if ( cs_gui_strcmp(model, "Rij-epsilon")
2894 || cs_gui_strcmp(model, "Rij-SSG")) {
2895
2896 cs_field_t *c_rij = cs_field_by_name_try("rij");
2897 cs_field_t *c_eps = cs_field_by_name("epsilon");
2898
2899 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
2900 cs_lnum_t c_id = cell_ids[e_id];
2901 for (int drij = 0; drij < 6; drij++) {
2902 c_rij->val[6*c_id + drij] = ini_vals[7*e_id + drij];
2903 c_eps->val[c_id] = ini_vals[7 * e_id + 6];
2904 }
2905 }
2906 }
2907 else if (cs_gui_strcmp(model, "Rij-EBRSM")) {
2908 cs_field_t *c_rij = cs_field_by_name_try("rij");
2909 cs_field_t *c_eps = cs_field_by_name("epsilon");
2910 cs_field_t *c_alp = cs_field_by_name("alpha");
2911
2912 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
2913 cs_lnum_t c_id = cell_ids[e_id];
2914 for (int drij = 0; drij < 6; drij++) {
2915 c_rij->val[6*c_id + drij] = ini_vals[8*e_id + drij];
2916 c_eps->val[c_id] = ini_vals[8 * e_id + 6];
2917 c_alp->val[c_id] = ini_vals[8 * e_id + 7];
2918 }
2919 }
2920 }
2921 else if (cs_gui_strcmp(model, "v2f-BL-v2/k")) {
2922
2923 cs_field_t *c_k = cs_field_by_name("k");
2924 cs_field_t *c_eps = cs_field_by_name("epsilon");
2925 cs_field_t *c_phi = cs_field_by_name("phi");
2926 cs_field_t *c_alp = cs_field_by_name("alpha");
2927
2928 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
2929 cs_lnum_t c_id = cell_ids[e_id];
2930
2931 c_k->val[c_id] = ini_vals[4 * e_id];
2932 c_eps->val[c_id] = ini_vals[4 * e_id + 1];
2933 c_phi->val[c_id] = ini_vals[4 * e_id + 2];
2934 c_alp->val[c_id] = ini_vals[4 * e_id + 3];
2935 }
2936 }
2937 else if (cs_gui_strcmp(model, "k-omega-SST")) {
2938
2939 cs_field_t *c_k = cs_field_by_name("k");
2940 cs_field_t *c_ome = cs_field_by_name("omega");
2941
2942 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
2943 cs_lnum_t c_id = cell_ids[e_id];
2944
2945 c_k->val[c_id] = ini_vals[2 * e_id];
2946 c_ome->val[c_id] = ini_vals[2 * e_id + 1];
2947 }
2948 }
2949 else if (cs_gui_strcmp(model, "Spalart-Allmaras")) {
2950 cs_field_t *c_nu = cs_field_by_name("nu_tilda");
2951
2952 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
2953 cs_lnum_t c_id = cell_ids[e_id];
2954 c_nu->val[c_id] = ini_vals[e_id];
2955 }
2956 }
2957
2958 else
2959 bft_error(__FILE__, __LINE__, 0,
2960 _("Invalid turbulence model: %s.\n"), model);
2961
2962 BFT_FREE(ini_vals);
2963
2964 }
2965 }
2966 }
2967
2968 /* Thermal scalar initialization */
2969 if (cs_gui_thermal_model() > 0) {
2970
2971 const char *formula_sca = NULL;
2972 cs_tree_node_t *tn_sca
2973 = cs_tree_get_node(cs_glob_tree,
2974 "thermophysical_models/"
2975 "thermal_scalar/variable/formula");
2976 tn_sca = _add_zone_id_test_attribute(tn_sca, z->id);
2977 formula_sca = cs_tree_node_get_value_str(tn_sca);
2978
2979 /* For non-specific physics defined with the GUI,
2980 the thermal variable can only be temperature or enthalpy
2981 (as the thermal model is on) */
2982
2983 cs_field_t *c = cs_thermal_model_field();
2984
2985 assert(c != NULL);
2986
2987 if (formula_sca != NULL) {
2988 cs_real_t *ini_vals = cs_meg_initialization(z, "thermal");
2989 if (ini_vals != NULL) {
2990 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
2991 cs_lnum_t c_id = cell_ids[e_id];
2992 c->val[c_id] = ini_vals[e_id];
2993 }
2994 BFT_FREE(ini_vals);
2995 }
2996 }
2997 /* If no formula was provided, the previous field values are
2998 kept (allowing mode-specific automatic initialization). */
2999 }
3000
3001 /* User Scalars initialization */
3002 int n_fields = cs_field_n_fields();
3003
3004 for (int f_id = 0; f_id < n_fields; f_id++) {
3005
3006 const cs_field_t *f = cs_field_by_id(f_id);
3007
3008 if ( f->type & CS_FIELD_USER
3009 && f->location_id == CS_MESH_LOCATION_CELLS) {
3010
3011 const char *formula_sca = NULL;
3012
3013 cs_tree_node_t *tn_sca = NULL;
3014 tn_sca = cs_tree_get_node(cs_glob_tree,
3015 "additional_scalars/variable");
3016 tn_sca = cs_tree_node_get_sibling_with_tag(tn_sca, "name", f->name);
3017 tn_sca = cs_tree_get_node(tn_sca, "formula");
3018 tn_sca = _add_zone_id_test_attribute(tn_sca, z->id);
3019 formula_sca = cs_tree_node_get_value_str(tn_sca);
3020
3021 if (formula_sca != NULL) {
3022 cs_real_t *ini_vals = cs_meg_initialization(z, f->name);
3023 if (ini_vals != NULL) {
3024 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
3025 cs_lnum_t c_id = cell_ids[e_id];
3026 f->val[c_id] = ini_vals[e_id];
3027 }
3028 BFT_FREE(ini_vals);
3029 }
3030 }
3031 }
3032 }
3033
3034 /* Meteo Scalars initialization */
3035 if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] > -1) {
3036
3037 cs_tree_node_t *tn_m0
3038 = cs_tree_get_node(cs_glob_tree,
3039 "thermophysical_models/atmospheric_flows");
3040
3041 const char *name = NULL;
3042 const char *formula_meteo = NULL;
3043
3044 int size = cs_tree_get_sub_node_count_simple(tn_m0, "variable");
3045
3046 for (int j = 0; j < size; j++) {
3047 cs_tree_node_t *tn_meteo = cs_tree_get_node(tn_m0, "variable");
3048 for (int i = 1;
3049 tn_meteo != NULL && i < j + 1;
3050 i++) {
3051 tn_meteo = cs_tree_node_get_next_of_name(tn_meteo);
3052 }
3053 cs_tree_node_t *tn_meteo2 = tn_meteo;
3054 tn_meteo = cs_tree_get_node(tn_meteo, "name");
3055 name = cs_tree_node_get_value_str(tn_meteo);
3056
3057 cs_field_t *c = cs_field_by_name_try(name);
3058
3059 snprintf(z_id_str, 31, "%d", z_id);
3060 const char *zone_id
3061 = cs_tree_node_get_child_value_str(tn_meteo2, "zone_id");
3062
3063 if (cs_gui_strcmp(zone_id, z_id_str))
3064 tn_meteo2 = cs_tree_get_node(tn_meteo2, "formula");
3065 else
3066 tn_meteo2 = NULL;
3067
3068 formula_meteo = cs_tree_node_get_value_str(tn_meteo2);
3069
3070 if (formula_meteo != NULL) {
3071 cs_real_t *ini_vals = cs_meg_initialization(z, c->name);
3072 if (ini_vals != NULL) {
3073 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
3074 cs_lnum_t c_id = cell_ids[e_id];
3075 c->val[c_id] = ini_vals[e_id];
3076 }
3077 BFT_FREE(ini_vals);
3078 }
3079 }
3080 else {
3081 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
3082 cs_lnum_t c_id = cell_ids[e_id];
3083 c->val[c_id] = 0.0;
3084 }
3085 }
3086 }
3087
3088 }
3089
3090 /* Combustion Scalars initialization */
3091 if (cs_glob_physical_model_flag[CS_COMBUSTION_3PT] > -1) {
3092
3093 cs_tree_node_t *tn_gas
3094 = cs_tree_get_node(cs_glob_tree,
3095 "thermophysical_models/gas_combustion");
3096
3097 const char *name = NULL;
3098 const char *formula_comb = NULL;
3099
3100 int size = cs_tree_get_sub_node_count_simple(tn_gas, "variable");
3101
3102 for (int j = 0; j < size; j++) {
3103 cs_tree_node_t *tn_combustion = cs_tree_get_node(tn_gas,
3104 "variable");
3105 for (int i = 1;
3106 tn_combustion != NULL && i < j + 1;
3107 i++) {
3108 tn_combustion = cs_tree_node_get_next_of_name(tn_combustion);
3109 }
3110
3111 cs_tree_node_t *tn_combustion2 = tn_combustion;
3112 tn_combustion = cs_tree_get_node(tn_combustion, "name");
3113 name = cs_tree_node_get_value_str(tn_combustion);
3114
3115 tn_combustion2 = cs_tree_get_node(tn_combustion2, "formula");
3116 tn_combustion2 = _add_zone_id_test_attribute(tn_combustion2, z->id);
3117
3118 cs_field_t *c_comb = cs_field_by_name_try(name);
3119
3120 formula_comb = cs_tree_node_get_value_str(tn_combustion2);
3121
3122 if (formula_comb != NULL) {
3123 cs_real_t *ini_vals = cs_meg_initialization(z, c_comb->name);
3124 if (ini_vals != NULL) {
3125 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
3126 cs_lnum_t c_id = cell_ids[e_id];
3127 c_comb->val[c_id] = ini_vals[e_id];
3128 }
3129 BFT_FREE(ini_vals);
3130 }
3131 }
3132 }
3133
3134 }
3135
3136 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
3137 const char *formula = NULL;
3138 const char *buff = NULL;
3139 const char *name[] = {"pressure", "temperature", "total_energy",
3140 "density"};
3141
3142 ccfth = 10000;
3143 for (int j = 0; j < 4; j++) {
3144
3145 cs_tree_node_t *tn = NULL;
3146 if (j < 3) {
3147 tn = cs_tree_find_node(cs_glob_tree, "variable");
3148 while (tn != NULL) {
3149 const char *name_tn
3150 = cs_tree_node_get_child_value_str(tn, "name");
3151 if (cs_gui_strcmp(name_tn, name[j]))
3152 break;
3153 else
3154 tn = cs_tree_find_node_next(cs_glob_tree, tn, "variable");
3155 }
3156 }
3157 else {
3158 tn = cs_tree_find_node(cs_glob_tree, "property");
3159 while (tn != NULL) {
3160 const char *name_tn
3161 = cs_tree_node_get_child_value_str(tn, "name");
3162 if (cs_gui_strcmp(name_tn, name[j]))
3163 break;
3164 else
3165 tn = cs_tree_find_node_next(cs_glob_tree, tn, "property");
3166 }
3167 }
3168 tn = cs_tree_get_node(tn, "formula");
3169 tn =_add_zone_id_test_attribute(tn, z->id);
3170 buff = cs_tree_node_get_child_value_str(tn, "status");
3171
3172 if (cs_gui_strcmp(buff, "on")) {
3173 if (j == 0)
3174 ccfth = ccfth * 2;
3175 else if (j == 1)
3176 ccfth = ccfth * 5;
3177 else if (j == 2)
3178 ccfth = ccfth * 7;
3179 else if (j == 3)
3180 ccfth = ccfth * 3;
3181
3182 cs_field_t *c = cs_field_by_name_try(name[j]);
3183
3184 formula = cs_tree_node_get_value_str(tn);
3185
3186 if (formula != NULL) {
3187 cs_real_t *ini_vals = cs_meg_initialization(z, c->name);
3188 if (ini_vals != NULL) {
3189 for (cs_lnum_t e_id = 0; e_id < n_cells; e_id++) {
3190 cs_lnum_t c_id = cell_ids[e_id];
3191 c->val[c_id] = ini_vals[e_id];
3192 }
3193 BFT_FREE(ini_vals);
3194 }
3195 }
3196 }
3197
3198 }
3199 *iccfth = ccfth;
3200 }
3201
3202 } /* END OF ISUITE == 0 */
3203 }
3204 } /* zones+1 */
3205 }
3206
3207 /*----------------------------------------------------------------------------
3208 * User law for material properties
3209 *
3210 * Fortran Interface:
3211 *
3212 * subroutine uiphyv
3213 * *****************
3214 *
3215 * integer iviscv <-- pointer for volumic viscosity viscv
3216 *----------------------------------------------------------------------------*/
3217
CS_PROCF(uiphyv,UIPHYV)3218 void CS_PROCF(uiphyv, UIPHYV)(const int *iviscv)
3219 {
3220 double time0 = cs_timer_wtime();
3221
3222 int n_zones_pp
3223 = cs_volume_zone_n_type_zones(CS_VOLUME_ZONE_PHYSICAL_PROPERTIES);
3224 int n_zones = cs_volume_zone_n_zones();
3225 /* law for density (built-in for all current integrated physical models) */
3226 if (cs_glob_fluid_properties->irovar == 1) {
3227 cs_field_t *c_rho = CS_F_(rho);
3228 if (n_zones_pp > 0) {
3229 for (int z_id = 0; z_id < n_zones; z_id++) {
3230 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
3231 if (z->type & CS_VOLUME_ZONE_PHYSICAL_PROPERTIES)
3232 _physical_property(c_rho, z);
3233 }
3234 }
3235 }
3236
3237 /* law for molecular viscosity */
3238 if (cs_glob_fluid_properties->ivivar == 1) {
3239 cs_field_t *c_mu = CS_F_(mu);
3240 if (n_zones_pp > 0) {
3241 for (int z_id = 0; z_id < n_zones; z_id++) {
3242 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
3243 if (z->type & CS_VOLUME_ZONE_PHYSICAL_PROPERTIES)
3244 _physical_property(c_mu, z);
3245 }
3246 }
3247 }
3248
3249 /* law for specific heat */
3250 if (cs_glob_fluid_properties->icp > 0) {
3251 cs_field_t *c_cp = CS_F_(cp);
3252 if (n_zones_pp > 0) {
3253 for (int z_id = 0; z_id < n_zones; z_id++) {
3254 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
3255 if (z->type & CS_VOLUME_ZONE_PHYSICAL_PROPERTIES)
3256 _physical_property(c_cp, z);
3257 }
3258 }
3259 }
3260
3261 /* law for thermal conductivity */
3262 if (cs_glob_thermal_model->itherm != CS_THERMAL_MODEL_NONE) {
3263
3264 cs_field_t *cond_dif = NULL;
3265
3266 cs_field_t *_th_f[] = {CS_F_(t), CS_F_(h), CS_F_(e_tot)};
3267
3268 for (int i = 0; i < 3; i++)
3269 if (_th_f[i]) {
3270 if ((_th_f[i])->type & CS_FIELD_VARIABLE) {
3271 int k = cs_field_key_id("diffusivity_id");
3272 int cond_diff_id = cs_field_get_key_int(_th_f[i], k);
3273 if (cond_diff_id > -1) {
3274 cond_dif = cs_field_by_id(cond_diff_id);
3275 if (n_zones_pp > 0) {
3276 for (int z_id = 0; z_id < n_zones; z_id++) {
3277 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
3278 if (z->type & CS_VOLUME_ZONE_PHYSICAL_PROPERTIES)
3279 _physical_property(cond_dif, z);
3280 }
3281 }
3282 }
3283 break;
3284 }
3285 }
3286 }
3287
3288 /* law for volumic viscosity (compressible model) */
3289 if (cs_glob_physical_model_flag[CS_COMPRESSIBLE] > -1) {
3290 if (*iviscv > 0) {
3291 cs_field_t *c = cs_field_by_name_try("volume_viscosity");
3292 if (n_zones_pp > 0 && c != NULL) {
3293 for (int z_id = 0; z_id < n_zones; z_id++) {
3294 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
3295 if (z->type & CS_VOLUME_ZONE_PHYSICAL_PROPERTIES)
3296 _physical_property(c, z);
3297 }
3298 }
3299 }
3300 }
3301
3302 /* law for scalar diffusivity */
3303 int user_id = -1;
3304 int n_fields = cs_field_n_fields();
3305 const int kivisl = cs_field_key_id("diffusivity_id");
3306 const int kscavr = cs_field_key_id("first_moment_id");
3307
3308 for (int f_id = 0; f_id < n_fields; f_id++) {
3309
3310 const cs_field_t *f = cs_field_by_id(f_id);
3311
3312 if ( (f->type & CS_FIELD_VARIABLE)
3313 && (f->type & CS_FIELD_USER)) {
3314 user_id++;
3315
3316 if ( cs_field_get_key_int(f, kscavr) < 0
3317 && cs_field_get_key_int(f, kivisl) >= 0) {
3318
3319 /* Get diffusivity pointer.
3320 * diff_id is >= 0 since it is a part of the if test
3321 * above!
3322 */
3323 int diff_id = cs_field_get_key_int(f, kivisl);
3324 cs_field_t *c_prop = cs_field_by_id(diff_id);
3325
3326 if (n_zones_pp > 0) {
3327 for (int z_id = 0; z_id < n_zones; z_id++) {
3328 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
3329 if (z->type & CS_VOLUME_ZONE_PHYSICAL_PROPERTIES) {
3330 const char *law = _property_formula(c_prop->name, z->name);
3331 if (law != NULL) {
3332 _physical_property(c_prop, z);
3333 if (cs_glob_fluid_properties->irovar == 1) {
3334 cs_real_t *c_rho = CS_F_(rho)->val;
3335 for (cs_lnum_t e_id = 0; e_id < z->n_elts; e_id++) {
3336 cs_lnum_t c_id = z->elt_ids[e_id];
3337 c_prop->val[c_id] *= c_rho[c_id];
3338 }
3339 }
3340 else {
3341 for (cs_lnum_t e_id = 0; e_id < z->n_elts; e_id++) {
3342 cs_lnum_t c_id = z->elt_ids[e_id];
3343 c_prop->val[c_id] *= cs_glob_fluid_properties->ro0;
3344 }
3345 }
3346 cs_gui_add_mei_time(cs_timer_wtime() - time0);
3347 #if _XML_DEBUG
3348 bft_printf("==> %s\n", __func__);
3349 bft_printf("--law for the diffusivity coefficient "
3350 "of the scalar '%s' over zone '%s':\n %s\n",
3351 f->name, z->name, law);
3352 #endif
3353 }
3354 }
3355 }
3356 }
3357 }
3358 }
3359 }
3360
3361 }
3362
3363 /*----------------------------------------------------------------------------
3364 * extra operations
3365 *
3366 * Fortran Interface:
3367 *
3368 * SUBROUTINE UIEXOP
3369 * *****************
3370 *
3371 *----------------------------------------------------------------------------*/
3372
CS_PROCF(uiexop,UIEXOP)3373 void CS_PROCF (uiexop, UIEXOP)(void)
3374 {
3375 cs_gui_balance_by_zone();
3376 cs_gui_pressure_drop_by_zone();
3377 }
3378
3379 /*----------------------------------------------------------------------------
3380 * groundwater model : read laws for capacity, saturation and permeability
3381 *
3382 * Fortran Interface:
3383 *
3384 * subroutine uidapp
3385 * *****************
3386 * integer permeability <-- permeability type
3387 * integer diffusion <-- diffusion type
3388 * integer gravity <-- check if gravity is taken into account
3389 * double gravity_x <-- gravity direction
3390 * double gravity_y <-- gravity direction
3391 * double gravity_z <-- gravity direction
3392 * integer unsaturated <-- unsaturated zone taken into account
3393 *----------------------------------------------------------------------------*/
3394
CS_PROCF(uidapp,UIDAPP)3395 void CS_PROCF (uidapp, UIDAPP) (const int *permeability,
3396 const int *diffusion,
3397 const int *gravity,
3398 const cs_real_t *gravity_x,
3399 const cs_real_t *gravity_y,
3400 const cs_real_t *gravity_z,
3401 const int *unsaturated)
3402 {
3403 const cs_real_3_t *restrict cell_cen
3404 = (const cs_real_3_t *restrict)cs_glob_mesh_quantities->cell_cen;
3405
3406 const cs_real_3_t *vel = (const cs_real_3_t *)(CS_F_(vel)->val);
3407
3408 cs_field_t *fsaturation = cs_field_by_name_try("saturation");
3409 cs_field_t *fcapacity = cs_field_by_name_try("capacity");
3410 cs_field_t *fpermeability = cs_field_by_name_try("permeability");
3411 cs_field_t *fhhead = CS_F_(head);
3412 cs_field_t *fsoil_density = cs_field_by_name_try("soil_density");
3413
3414 cs_real_t *saturation_field = fsaturation->val;
3415 cs_real_t *capacity_field = fcapacity->val;
3416 cs_real_t *h_head_field = fhhead->val;
3417 cs_real_t *soil_density = fsoil_density->val;
3418
3419 cs_real_t *permeability_field = NULL;
3420 cs_real_6_t *permeability_field_v = NULL;
3421
3422 cs_gnum_t cw[3];
3423
3424 if (*permeability == 0)
3425 permeability_field = fpermeability->val;
3426 else
3427 permeability_field_v = (cs_real_6_t *)fpermeability->val;
3428
3429 cs_tree_node_t *tn_gw
3430 = cs_tree_get_node(cs_glob_tree, "thermophysical_models/groundwater");
3431
3432 /* number of volumic zone */
3433
3434 int n_zones = cs_volume_zone_n_zones();
3435
3436 for (int z_id = 0; z_id < n_zones; z_id++) {
3437
3438 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
3439
3440 if (_zone_id_is_type(z->id, "groundwater_law")) {
3441
3442 cs_tree_node_t *tn_zl = cs_tree_get_node(tn_gw, "groundwater_law");
3443 tn_zl = _add_zone_id_test_attribute(tn_zl, z->id);
3444
3445 const cs_lnum_t n_cells = z->n_elts;
3446 const cs_lnum_t *cell_ids = z->elt_ids;
3447
3448 char z_id_str[32];
3449 snprintf(z_id_str, 31, "%d", z_id);
3450
3451 /* get ground properties for each zone */
3452
3453 /* get soil density by zone */
3454 cs_real_t rhosoil = 0.;
3455
3456 cs_gui_node_get_child_real(tn_zl, "soil_density", &rhosoil);
3457
3458 for (cs_lnum_t icel = 0; icel < n_cells; icel++) {
3459 cs_lnum_t iel = cell_ids[icel];
3460 soil_density[iel] = rhosoil;
3461 }
3462
3463 const char *mdl = cs_tree_node_get_child_value_str(tn_zl, "model");
3464
3465 /* law for permeability */
3466 /* TODO: rename it in GUI, it is not Van Genuchten if saturated */
3467
3468 if (cs_gui_strcmp(mdl, "VanGenuchten")) {
3469
3470 cs_real_t alpha_param, ks_param, l_param, n_param;
3471 cs_real_t thetas_param, thetar_param;
3472 cs_real_t ks_xx, ks_yy, ks_zz, ks_xy, ks_xz, ks_yz;
3473
3474 cs_tree_node_t *tn_m
3475 = cs_tree_node_get_child(tn_zl, "VanGenuchten_parameters");
3476
3477 /* Van Genuchten parameters */
3478 if (*unsaturated) {
3479 cs_gui_node_get_child_real(tn_m, "alpha", &alpha_param);
3480 cs_gui_node_get_child_real(tn_m, "l", &l_param);
3481 cs_gui_node_get_child_real(tn_m, "n", &n_param);
3482 cs_gui_node_get_child_real(tn_m, "thetar", &thetar_param);
3483 }
3484
3485 cs_gui_node_get_child_real(tn_m, "thetas", &thetas_param);
3486
3487 if (*permeability == 0)
3488 cs_gui_node_get_child_real(tn_m, "ks", &ks_param);
3489 else {
3490 cs_gui_node_get_child_real(tn_m, "ks_xx", &ks_xx);
3491 cs_gui_node_get_child_real(tn_m, "ks_yy", &ks_yy);
3492 cs_gui_node_get_child_real(tn_m, "ks_zz", &ks_zz);
3493 cs_gui_node_get_child_real(tn_m, "ks_xy", &ks_xy);
3494 cs_gui_node_get_child_real(tn_m, "ks_yz", &ks_yz);
3495 cs_gui_node_get_child_real(tn_m, "ks_xz", &ks_xz);
3496 }
3497
3498 /* unsaturated zone considered */
3499 if (*unsaturated) {
3500 for (cs_lnum_t icel = 0; icel < n_cells; icel++) {
3501 cs_lnum_t iel = cell_ids[icel];
3502 cs_real_t head = h_head_field[iel];
3503
3504 if (*gravity == 1)
3505 head -= ( cell_cen[iel][0] * *gravity_x
3506 + cell_cen[iel][1] * *gravity_y
3507 + cell_cen[iel][2] * *gravity_z);
3508
3509 if (head >= 0) {
3510 capacity_field[iel] = 0.;
3511 saturation_field[iel] = thetas_param;
3512
3513 if (*permeability == 0)
3514 permeability_field[iel] = ks_param;
3515 else {
3516 permeability_field_v[iel][0] = ks_xx;
3517 permeability_field_v[iel][1] = ks_yy;
3518 permeability_field_v[iel][2] = ks_zz;
3519 permeability_field_v[iel][3] = ks_xy;
3520 permeability_field_v[iel][4] = ks_yz;
3521 permeability_field_v[iel][5] = ks_xz;
3522 }
3523 }
3524 else {
3525 cs_real_t m_param = 1 - 1 / n_param;
3526 cs_real_t tmp1 = pow(fabs(alpha_param * head), n_param);
3527 cs_real_t tmp2 = 1. / (1. + tmp1);
3528 cs_real_t se_param = pow(tmp2, m_param);
3529 cs_real_t perm = pow(se_param, l_param) *
3530 pow((1. - pow((1. - tmp2), m_param)), 2);
3531
3532 capacity_field[iel] = -m_param * n_param * tmp1 *
3533 (thetas_param - thetar_param) *
3534 se_param * tmp2 / head;
3535 saturation_field[iel] = thetar_param +
3536 se_param * (thetas_param - thetar_param);
3537
3538 if (*permeability == 0)
3539 permeability_field[iel] = perm * ks_param;
3540 else {
3541 permeability_field_v[iel][0] = perm * ks_xx;
3542 permeability_field_v[iel][1] = perm * ks_yy;
3543 permeability_field_v[iel][2] = perm * ks_zz;
3544 permeability_field_v[iel][3] = perm * ks_xy;
3545 permeability_field_v[iel][4] = perm * ks_yz;
3546 permeability_field_v[iel][5] = perm * ks_xz;
3547 }
3548 }
3549 }
3550 }
3551 else { /* saturated */
3552 for (cs_lnum_t icel = 0; icel < n_cells; icel++) {
3553 cs_lnum_t iel = cell_ids[icel];
3554 capacity_field[iel] = 0.;
3555 saturation_field[iel] = thetas_param;
3556
3557 if (*permeability == 0)
3558 permeability_field[iel] = ks_param;
3559 else {
3560 permeability_field_v[iel][0] = ks_xx;
3561 permeability_field_v[iel][1] = ks_yy;
3562 permeability_field_v[iel][2] = ks_zz;
3563 permeability_field_v[iel][3] = ks_xy;
3564 permeability_field_v[iel][4] = ks_yz;
3565 permeability_field_v[iel][5] = ks_xz;
3566 }
3567 }
3568 }
3569
3570 } else {
3571 /* user law for permeability */
3572 const char *formula
3573 = cs_tree_node_get_child_value_str(tn_zl, "formula");
3574
3575 if (formula != NULL) {
3576 cs_field_t *fmeg[3] = {fcapacity, fsaturation, fpermeability};
3577 cs_meg_volume_function(z, fmeg);
3578 }
3579 }
3580
3581 const int kivisl = cs_field_key_id("diffusivity_id");
3582 int n_fields = cs_field_n_fields();
3583
3584
3585 /* get diffusivity and Kd for each scalar defined by the user on
3586 current zone (and kplus and kminus only for scalars
3587 with kinetic model) */
3588 for (int f_id = 0; f_id < n_fields; f_id++) {
3589
3590 cs_field_t *f = cs_field_by_id(f_id);
3591
3592 if ( (f->type & CS_FIELD_VARIABLE)
3593 && (f->type & CS_FIELD_USER)) {
3594
3595 /* get kd for current scalar and current zone */
3596 char *kdname = NULL;
3597 int len = strlen(f->name) + 4;
3598 BFT_MALLOC(kdname, len, char);
3599 strcpy(kdname, f->name);
3600 strcat(kdname, "_kd");
3601 cs_field_t *fkd = cs_field_by_name_try(kdname);
3602 BFT_FREE(kdname);
3603
3604 cs_real_t kd_val = 0., diff_val = 0.;
3605 cs_tree_node_t *tn_s = cs_tree_get_node(tn_zl, "scalar");
3606 tn_s = cs_tree_node_get_sibling_with_tag(tn_s, "name", f->name);
3607
3608 cs_gui_node_get_child_real(tn_s, "kd", &kd_val);
3609 cs_gui_node_get_child_real(tn_s, "diffusivity", &diff_val);
3610
3611 for (cs_lnum_t icel = 0; icel < n_cells; icel++) {
3612 cs_lnum_t iel = cell_ids[icel];
3613 fkd->val[iel] = kd_val;
3614 }
3615
3616 /* get diffusivity for current scalar and current zone */
3617 int diff_id = cs_field_get_key_int(f, kivisl);
3618 cs_field_t *fdiff = cs_field_by_id(diff_id);
3619
3620 for (cs_lnum_t icel = 0; icel < n_cells; icel++) {
3621 cs_lnum_t iel = cell_ids[icel];
3622 fdiff->val[iel] = saturation_field[iel]*diff_val;
3623 }
3624
3625 /* get kplus and kminus for current scalar and current zone
3626 (if EK model is chosen) */
3627 cs_gwf_soilwater_partition_t sorption_scal;
3628 int key_part = cs_field_key_id("gwf_soilwater_partition");
3629 cs_field_t *kp, *km;
3630 cs_field_get_key_struct(f, key_part, &sorption_scal);
3631
3632 if (sorption_scal.kinetic == 1) {
3633
3634 kp = cs_field_by_id(sorption_scal.ikp);
3635 km = cs_field_by_id(sorption_scal.ikm);
3636
3637 cs_real_t kp_val = 0., km_val = 0.;
3638 cs_gui_node_get_child_real(tn_s, "kplus", &kp_val);
3639 cs_gui_node_get_child_real(tn_s, "kminus", &km_val);
3640
3641 for (cs_lnum_t icel = 0; icel < n_cells; icel++) {
3642 cs_lnum_t iel = cell_ids[icel];
3643 kp->val[iel] = kp_val;
3644 km->val[iel] = km_val;
3645 }
3646
3647 }
3648
3649 }
3650 }
3651
3652 /* get dispersion coefficient */
3653 if (*diffusion == 1) { /* anisotropic dispersion */
3654 /* TODO use a dedicated tensor field by species */
3655
3656 cs_field_t *fturbvisco
3657 = cs_field_by_name_try("anisotropic_turbulent_viscosity");
3658 cs_real_6_t *visten_v = (cs_real_6_t *)fturbvisco->val;
3659
3660 cs_real_t long_diffus;
3661 double trans_diffus;
3662
3663 cs_tree_node_t *tn
3664 = cs_tree_node_get_child(tn_zl, "diffusion_coefficient");
3665 cs_gui_node_get_child_real(tn, "longitudinal", &long_diffus);
3666 cs_gui_node_get_child_real(tn, "transverse", &trans_diffus);
3667
3668 for (cs_lnum_t icel = 0; icel < n_cells; icel++) {
3669 cs_lnum_t iel = cell_ids[icel];
3670 double norm = sqrt(vel[iel][0] * vel[iel][0] +
3671 vel[iel][1] * vel[iel][1] +
3672 vel[iel][2] * vel[iel][2]);
3673 double tmp = trans_diffus * norm;
3674 double diff = long_diffus - trans_diffus;
3675 double denom = norm + 1.e-15;
3676 visten_v[iel][0] = tmp + diff * vel[iel][0] * vel[iel][0] / denom;
3677 visten_v[iel][1] = tmp + diff * vel[iel][1] * vel[iel][1] / denom;
3678 visten_v[iel][2] = tmp + diff * vel[iel][2] * vel[iel][2] / denom;
3679 visten_v[iel][3] = diff * vel[iel][1] * vel[iel][0] / denom;
3680 visten_v[iel][4] = diff * vel[iel][1] * vel[iel][2] / denom;
3681 visten_v[iel][5] = diff * vel[iel][2] * vel[iel][0] / denom;
3682 }
3683
3684 }
3685 else { /* isotropic dispersion */
3686 /* - same value of isotropic dispersion for each species
3687 - assigned to diffusivity field of each species
3688 TODO: allow to specifiy one value by species in GUI */
3689
3690 double diffus;
3691
3692 cs_tree_node_t *tn = cs_tree_node_get_child(tn_zl,
3693 "diffusion_coefficient");
3694 cs_gui_node_get_child_real(tn, "isotropic", &diffus);
3695
3696 for (int f_id = 0; f_id < n_fields; f_id++) {
3697 cs_field_t *f = cs_field_by_id(f_id);
3698 if ( (f->type & CS_FIELD_VARIABLE)
3699 && (f->type & CS_FIELD_USER)) {
3700 int diff_id = cs_field_get_key_int(f, kivisl);
3701 cs_field_t *fdiff = cs_field_by_id(diff_id);
3702 cs_real_t *visten = fdiff->val;
3703
3704 /* WARNING: dispersion adds up to diffusivity
3705 already assigned above */
3706 for (cs_lnum_t l_id = 0; l_id < n_cells; l_id++) {
3707 cs_lnum_t iel = cell_ids[l_id];
3708 cs_real_t norm = sqrt(vel[iel][0] * vel[iel][0] +
3709 vel[iel][1] * vel[iel][1] +
3710 vel[iel][2] * vel[iel][2]);
3711 visten[iel] = visten[iel] + diffus * norm;
3712 }
3713 }
3714 }
3715
3716 }
3717 }
3718 }
3719
3720 /* check values */
3721
3722 {
3723 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
3724
3725 cw[0] = 0; cw[1] = 0; cw[2] = 0;
3726
3727 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
3728 if (saturation_field[iel] > 1. || saturation_field[iel] < 0.)
3729 cw[0] += 1;
3730
3731 if (capacity_field[iel] < 0.)
3732 cw[1] += 1;
3733
3734 if (*permeability == 0) {
3735 if (permeability_field[iel] < 0.)
3736 cw[2] += 1;
3737 }
3738 }
3739
3740 cs_parall_counter(cw, 3);
3741
3742 if (cw[0] > 0)
3743 bft_printf(_("soil_tracer_law, WARNING:\n"
3744 " saturation is outside [0, 1] in %llu cells.\n"),
3745 (unsigned long long)(cw[0]));
3746
3747 if (cw[1] > 0)
3748 bft_printf(_("soil_tracer_law, WARNING:\n"
3749 " capacity is < 0 in %llu cells.\n"),
3750 (unsigned long long)(cw[1]));
3751
3752 if (cw[2] > 0)
3753 bft_printf(_("soil_tracer_law, WARNING:\n"
3754 " isotropic permeability is < 0 in %llu cells.\n"),
3755 (unsigned long long)(cw[2]));
3756 }
3757 }
3758
3759 /*----------------------------------------------------------------------------
3760 * Define fans with GUI
3761 *
3762 * Fortran Interface:
3763 *
3764 * SUBROUTINE UIFANS
3765 * *****************
3766 *
3767 *----------------------------------------------------------------------------*/
3768
CS_PROCF(uifans,UIFANS)3769 void CS_PROCF (uifans, UIFANS) (void)
3770 {
3771 cs_gui_define_fans();
3772 }
3773
3774 /*----------------------------------------------------------------------------
3775 * Define error estimators
3776 *
3777 * Fortran Interface:
3778 *
3779 * SUBROUTINE UIERES
3780 * *****************
3781 *
3782 *----------------------------------------------------------------------------*/
3783
CS_PROCF(uieres,UIERES)3784 void CS_PROCF (uieres, UIERES) (int *iescal,
3785 int *iespre,
3786 int *iesder,
3787 int *iescor,
3788 int *iestot)
3789 {
3790 cs_gui_error_estimator(iescal, iespre, iesder, iescor, iestot);
3791 }
3792
3793 /*============================================================================
3794 * Public function definitions
3795 *============================================================================*/
3796
3797 /*----------------------------------------------------------------------------
3798 * Initialize GUI reader structures.
3799 *----------------------------------------------------------------------------*/
3800
3801 /*-----------------------------------------------------------------------------
3802 * Free memory: clean global private variables.
3803 *----------------------------------------------------------------------------*/
3804
3805 void
cs_gui_finalize(void)3806 cs_gui_finalize(void)
3807 {
3808 cs_gui_boundary_conditions_free_memory();
3809 }
3810
3811 /*----------------------------------------------------------------------------*/
3812 /*!
3813 * \brief Compute GUI-defined head losses for a given volume zone.
3814 *
3815 * Head loss tensor coefficients for each cell are organized as follows:
3816 * cku11, cku22, cku33, cku12, cku13, cku23.
3817 *
3818 * \param[in] zone pointer to zone structure
3819 * \param[in] cvara_vel velocity values at the previous time step
3820 * \param[in, out] cku head loss coefficients
3821 */
3822 /*----------------------------------------------------------------------------*/
3823
3824 void
cs_gui_head_losses(const cs_zone_t * zone,const cs_real_3_t * cvara_vel,cs_real_t cku[][6])3825 cs_gui_head_losses(const cs_zone_t *zone,
3826 const cs_real_3_t *cvara_vel,
3827 cs_real_t cku[][6])
3828 {
3829 if (! (zone->type & CS_VOLUME_ZONE_HEAD_LOSS))
3830 return;
3831
3832 double c11, c12, c13, c21, c22, c23, c31, c32, c33;
3833
3834 const cs_lnum_t n_cells = zone->n_elts;
3835 const cs_lnum_t *cell_ids = zone->elt_ids;
3836
3837 char z_id_str[32];
3838 snprintf(z_id_str, 31, "%d", zone->id);
3839
3840 cs_tree_node_t *tn
3841 = cs_tree_get_node(cs_glob_tree,
3842 "thermophysical_models/head_losses/head_loss");
3843 tn = cs_tree_node_get_sibling_with_tag(tn, "zone_id", z_id_str);
3844
3845 double k11 = _c_head_losses(tn, "kxx");
3846 double k22 = _c_head_losses(tn, "kyy");
3847 double k33 = _c_head_losses(tn, "kzz");
3848
3849 double a11 = _c_head_losses(tn, "a11");
3850 double a12 = _c_head_losses(tn, "a12");
3851 double a13 = _c_head_losses(tn, "a13");
3852 double a21 = _c_head_losses(tn, "a21");
3853 double a22 = _c_head_losses(tn, "a22");
3854 double a23 = _c_head_losses(tn, "a23");
3855 double a31 = _c_head_losses(tn, "a31");
3856 double a32 = _c_head_losses(tn, "a32");
3857 double a33 = _c_head_losses(tn, "a33");
3858
3859 if ( cs_gui_is_equal_real(a12, 0.0)
3860 && cs_gui_is_equal_real(a13, 0.0)
3861 && cs_gui_is_equal_real(a23, 0.0)) {
3862
3863 c11 = k11;
3864 c22 = k22;
3865 c33 = k33;
3866 c12 = 0.0;
3867 c13 = 0.0;
3868 c23 = 0.0;
3869
3870 }
3871 else
3872 _matrix_base_conversion(a11, a12, a13, a21, a22, a23, a31, a32, a33,
3873 k11, 0.0, 0.0, 0.0, k22, 0.0, 0.0, 0.0, k33,
3874 &c11, &c12, &c13,
3875 &c21, &c22, &c23,
3876 &c31, &c32, &c33);
3877
3878 for (cs_lnum_t j = 0; j < n_cells; j++) {
3879 cs_lnum_t c_id = cell_ids[j];
3880 cs_real_t v = cs_math_3_norm(cvara_vel[c_id]);
3881 cku[j][0] = 0.5 * c11 * v;
3882 cku[j][1] = 0.5 * c22 * v;
3883 cku[j][2] = 0.5 * c33 * v;
3884 cku[j][3] = 0.5 * c12 * v;
3885 cku[j][4] = 0.5 * c23 * v;
3886 cku[j][5] = 0.5 * c13 * v;
3887 }
3888 }
3889
3890 /*-----------------------------------------------------------------------------
3891 * Selection of linear solvers.
3892 *----------------------------------------------------------------------------*/
3893
3894 void
cs_gui_linear_solvers(void)3895 cs_gui_linear_solvers(void)
3896 {
3897 bool multigrid = false;
3898 cs_sles_it_type_t sles_it_type = CS_SLES_N_IT_TYPES;
3899 cs_multigrid_type_t mg_type = CS_MULTIGRID_N_TYPES;
3900
3901 const char* algo_choice = NULL;
3902 const char* precond_choice = NULL;
3903
3904 const int n_max_iter_default = 10000;
3905
3906 int n_fields = cs_field_n_fields();
3907
3908 for (int f_id = 0; f_id < n_fields; f_id++) {
3909 cs_field_t *f = cs_field_by_id(f_id);
3910 if (f->type & CS_FIELD_VARIABLE) {
3911
3912 const char *ref_name = f->name;
3913
3914 if ( cs_gui_strcmp(f->name, "r11")
3915 || cs_gui_strcmp(f->name, "r22")
3916 || cs_gui_strcmp(f->name, "r33")
3917 || cs_gui_strcmp(f->name, "r12")
3918 || cs_gui_strcmp(f->name, "r23")
3919 || cs_gui_strcmp(f->name, "r13"))
3920 ref_name = "rij";
3921
3922 cs_tree_node_t *tn_v = _find_node_variable(ref_name);
3923
3924 int n_max_iter = n_max_iter_default;
3925 cs_gui_node_get_child_int(tn_v, "max_iter_number", &n_max_iter);
3926
3927 multigrid = false;
3928 sles_it_type = CS_SLES_N_IT_TYPES;
3929
3930 algo_choice = _variable_choice(tn_v, "solver_choice");
3931 precond_choice = _variable_choice(tn_v, "preconditioning_choice");
3932
3933 if (cs_gui_strcmp(algo_choice, "multigrid_k_cycle")) {
3934 multigrid = true;
3935 mg_type = CS_MULTIGRID_K_CYCLE;
3936 }
3937 else if (cs_gui_strcmp(algo_choice, "multigrid")) {
3938 multigrid = true;
3939 mg_type = CS_MULTIGRID_V_CYCLE;
3940 }
3941 else if (cs_gui_strcmp(algo_choice, "conjugate_gradient"))
3942 sles_it_type = CS_SLES_PCG;
3943 else if (cs_gui_strcmp(algo_choice, "flexible_conjugate_gradient"))
3944 sles_it_type = CS_SLES_FCG;
3945 else if (cs_gui_strcmp(algo_choice, "inexact_conjugate_gradient"))
3946 sles_it_type = CS_SLES_IPCG;
3947 else if (cs_gui_strcmp(algo_choice, "jacobi"))
3948 sles_it_type = CS_SLES_JACOBI;
3949 else if (cs_gui_strcmp(algo_choice, "bi_cgstab"))
3950 sles_it_type = CS_SLES_BICGSTAB;
3951 else if (cs_gui_strcmp(algo_choice, "bi_cgstab2"))
3952 sles_it_type = CS_SLES_BICGSTAB2;
3953 else if (cs_gui_strcmp(algo_choice, "gmres"))
3954 sles_it_type = CS_SLES_GMRES;
3955 else if (cs_gui_strcmp(algo_choice, "gauss_seidel"))
3956 sles_it_type = CS_SLES_P_GAUSS_SEIDEL;
3957 else if (cs_gui_strcmp(algo_choice, "symmetric_gauss_seidel"))
3958 sles_it_type = CS_SLES_P_SYM_GAUSS_SEIDEL;
3959 else if (cs_gui_strcmp(algo_choice, "PCR3"))
3960 sles_it_type = CS_SLES_PCR3;
3961
3962 /* If choice is "automatic" or unspecified, delay
3963 choice to cs_sles_default, so do nothing here */
3964
3965 if (sles_it_type < CS_SLES_N_IT_TYPES) {
3966
3967 int poly_degree = 0;
3968 bool pc_multigrid = false;
3969
3970 if (cs_gui_strcmp(precond_choice, "jacobi"))
3971 poly_degree = 0;
3972 else if (cs_gui_strcmp(precond_choice, "none"))
3973 poly_degree = -1;
3974 else if (cs_gui_strcmp(precond_choice, "polynomial"))
3975 poly_degree = 1;
3976 else if (cs_gui_strcmp(precond_choice, "multigrid_k_cycle")) {
3977 pc_multigrid = true;
3978 mg_type = CS_MULTIGRID_K_CYCLE;
3979 poly_degree = -1;
3980 }
3981 else if (cs_gui_strcmp(precond_choice, "multigrid_k_cycle_hpc")) {
3982 pc_multigrid = true;
3983 mg_type = CS_MULTIGRID_K_CYCLE_HPC;
3984 poly_degree = -1;
3985 }
3986 else if (cs_gui_strcmp(precond_choice, "multigrid")) {
3987 pc_multigrid = true;
3988 mg_type = CS_MULTIGRID_V_CYCLE;
3989 poly_degree = -1;
3990 }
3991 else { /* "automatic" or undefined */
3992 if (sles_it_type == CS_SLES_PCG) {
3993 pc_multigrid = true;
3994 mg_type = CS_MULTIGRID_V_CYCLE;
3995 poly_degree = -1;
3996 }
3997 }
3998
3999 cs_sles_it_t *c = cs_sles_it_define(f->id, NULL, sles_it_type,
4000 poly_degree, n_max_iter);
4001
4002 if (pc_multigrid) {
4003 cs_sles_pc_t *pc = cs_multigrid_pc_create(mg_type);
4004 cs_sles_it_transfer_pc(c, &pc);
4005 }
4006
4007 }
4008
4009 else if (multigrid == true) {
4010 cs_multigrid_t *mg = cs_multigrid_define(f->id, NULL, mg_type);
4011
4012 /* If we have convection, set appropriate options */
4013 if (f_id >= 0) {
4014 cs_var_cal_opt_t var_cal_opt;
4015 cs_field_get_key_struct(cs_field_by_id(f_id),
4016 cs_field_key_id("var_cal_opt"),
4017 &var_cal_opt);
4018 if (var_cal_opt.iconv > 0)
4019 cs_multigrid_set_solver_options
4020 (mg,
4021 CS_SLES_P_SYM_GAUSS_SEIDEL,
4022 CS_SLES_P_SYM_GAUSS_SEIDEL,
4023 CS_SLES_P_SYM_GAUSS_SEIDEL,
4024 100, /* n max cycles */
4025 3, /* n max iter for descent (default 2) */
4026 2, /* n max iter for ascent (default 10) */
4027 100, /* n max iter coarse solver */
4028 0, 0, 0, /* precond degree */
4029 -1, -1, 1); /* precision multiplier */
4030 }
4031 }
4032 }
4033 }
4034 }
4035
4036 /*-----------------------------------------------------------------------------
4037 * Define parallel IO settings.
4038 *----------------------------------------------------------------------------*/
4039
4040 void
cs_gui_parallel_io(void)4041 cs_gui_parallel_io(void)
4042 {
4043 int op_id;
4044 int rank_step = 0, block_size = -1;
4045
4046 cs_file_mode_t op_mode[2] = {CS_FILE_MODE_READ, CS_FILE_MODE_WRITE};
4047 const char *op_name[2] = {"read_method", "write_method"};
4048
4049 cs_tree_node_t *tn_bio = cs_tree_get_node(cs_glob_tree,
4050 "calculation_management/block_io");
4051 /* Block IO read and write method */
4052
4053 for (op_id = 0; op_id < 2; op_id++) {
4054
4055 cs_file_access_t m = CS_FILE_DEFAULT;
4056 const char *method_name
4057 = cs_tree_node_get_child_value_str(tn_bio, op_name[op_id]);
4058
4059 if (method_name != NULL) {
4060 if (!strcmp(method_name, "default"))
4061 m = CS_FILE_DEFAULT;
4062 else if (!strcmp(method_name, "stdio serial"))
4063 m = CS_FILE_STDIO_SERIAL;
4064 else if (!strcmp(method_name, "stdio parallel"))
4065 m = CS_FILE_STDIO_PARALLEL;
4066 else if (!strcmp(method_name, "mpi independent"))
4067 m = CS_FILE_MPI_INDEPENDENT;
4068 else if (!strcmp(method_name, "mpi noncollective"))
4069 m = CS_FILE_MPI_NON_COLLECTIVE;
4070 else if (!strcmp(method_name, "mpi collective"))
4071 m = CS_FILE_MPI_COLLECTIVE;
4072 #if defined(HAVE_MPI)
4073 cs_file_set_default_access(op_mode[op_id], m, MPI_INFO_NULL);
4074 #else
4075 cs_file_set_default_access(op_mode[op_id], m);
4076 #endif
4077 }
4078
4079 }
4080
4081 #if defined(HAVE_MPI)
4082
4083 /* Rank step and block buffer size */
4084
4085 cs_gui_node_get_child_int(tn_bio, "rank_step", &rank_step);
4086 cs_gui_node_get_child_int(tn_bio, "min_block_size", &block_size);
4087
4088 if (rank_step > 0 || block_size > -1) {
4089 int def_rank_step;
4090 cs_file_get_default_comm(&def_rank_step, NULL, NULL);
4091 size_t def_block_size = cs_parall_get_min_coll_buf_size();
4092 if (rank_step < 1)
4093 rank_step = def_rank_step;
4094 if (block_size < 0)
4095 block_size = def_block_size;
4096 cs_file_set_default_comm(rank_step, cs_glob_mpi_comm);
4097 cs_parall_set_min_coll_buf_size(block_size);
4098 }
4099
4100 #endif /* defined(HAVE_MPI) */
4101 }
4102
4103 /*-----------------------------------------------------------------------------
4104 * Set partitioning options.
4105 *----------------------------------------------------------------------------*/
4106
4107 void
cs_gui_partition(void)4108 cs_gui_partition(void)
4109 {
4110 cs_partition_algorithm_t a = CS_PARTITION_DEFAULT;
4111 bool ignore_perio = false;
4112 int rank_step = 1;
4113 int write_level = 1;
4114 int n_add_parts = 0;
4115 int *add_parts = NULL;
4116
4117 cs_tree_node_t *tn_p
4118 = cs_tree_get_node(cs_glob_tree, "calculation_management/partitioning");
4119
4120 /* Partitioning type */
4121 const char *part_name = cs_tree_node_get_child_value_str(tn_p, "type");
4122
4123 if (part_name != NULL) {
4124 if (!strcmp(part_name, "default"))
4125 a = CS_PARTITION_DEFAULT;
4126 else if (!strcmp(part_name, "morton sfc"))
4127 a = CS_PARTITION_SFC_MORTON_BOX;
4128 else if (!strcmp(part_name, "morton sfc cube"))
4129 a = CS_PARTITION_SFC_MORTON_CUBE;
4130 else if (!strcmp(part_name, "hilbert sfc"))
4131 a = CS_PARTITION_SFC_HILBERT_BOX;
4132 else if (!strcmp(part_name, "hilbert sfc cube"))
4133 a = CS_PARTITION_SFC_HILBERT_CUBE;
4134 else if (!strcmp(part_name, "scotch"))
4135 a = CS_PARTITION_SCOTCH;
4136 else if (!strcmp(part_name, "metis"))
4137 a = CS_PARTITION_METIS;
4138 else if (!strcmp(part_name, "block"))
4139 a = CS_PARTITION_BLOCK;
4140 }
4141
4142 /* Rank step */
4143
4144 cs_gui_node_get_child_int(tn_p, "rank_step", &rank_step);
4145
4146 /* Ignore periodicity option */
4147
4148 cs_gui_node_get_child_status_bool(tn_p, "ignore_periodicity", &ignore_perio);
4149
4150 /* Output option */
4151
4152 const char *s_output = cs_tree_node_get_child_value_str(tn_p, "output");
4153
4154 if (s_output != NULL) {
4155 if (!strcmp(s_output, "no"))
4156 write_level = 0;
4157 else if (!strcmp(s_output, "default"))
4158 write_level = 1;
4159 else if (!strcmp(s_output, "yes"))
4160 write_level = 2;
4161 }
4162
4163 /* List of partitions to output */
4164
4165 const char *s_list = cs_tree_node_get_child_value_str(tn_p, "partition_list");
4166
4167 if (s_list != NULL) {
4168 char *buf;
4169 BFT_MALLOC(buf, strlen(s_list)+1, char);
4170 strcpy(buf, s_list);
4171 char *p = strtok(buf, " \t,;");
4172 while (p != NULL) {
4173 int np = atoi(p);
4174 if (np > 1) {
4175 BFT_REALLOC(add_parts, n_add_parts + 1, int);
4176 add_parts[n_add_parts] = np;
4177 n_add_parts += 1;
4178 }
4179 p = strtok(NULL, " \t,;");
4180 }
4181 BFT_FREE(buf);
4182 }
4183
4184 /* Set options */
4185
4186 cs_partition_set_algorithm(CS_PARTITION_MAIN,
4187 a,
4188 rank_step,
4189 ignore_perio);
4190
4191 cs_partition_set_write_level(write_level);
4192
4193 if (n_add_parts > 0) {
4194 cs_partition_add_partitions(n_add_parts, add_parts);
4195 BFT_FREE(add_parts);
4196 }
4197 }
4198
4199 /*-----------------------------------------------------------------------------
4200 * Set MPI related algorithm options
4201 *----------------------------------------------------------------------------*/
4202
4203 void
cs_gui_mpi_algorithms(void)4204 cs_gui_mpi_algorithms(void)
4205 {
4206 cs_all_to_all_type_t a = CS_ALL_TO_ALL_MPI_DEFAULT;
4207
4208 cs_tree_node_t *tn_cm
4209 = cs_tree_get_node(cs_glob_tree, "calculation_management");
4210
4211 /* Partitioning type */
4212 const char *all_to_all_name
4213 = cs_tree_node_get_child_value_str(tn_cm, "all_to_all");
4214
4215 if (all_to_all_name != NULL) {
4216 if (!strcmp(all_to_all_name, "default"))
4217 a = CS_ALL_TO_ALL_MPI_DEFAULT;
4218 else if (!strcmp(all_to_all_name, "crystal router"))
4219 a = CS_ALL_TO_ALL_CRYSTAL_ROUTER;
4220 cs_all_to_all_set_type(a);
4221 }
4222 }
4223
4224 /*----------------------------------------------------------------------------
4225 * Determine porosity model type
4226 *----------------------------------------------------------------------------*/
4227
4228 void
cs_gui_porous_model(void)4229 cs_gui_porous_model(void)
4230 {
4231 int n_zones = cs_volume_zone_n_zones();
4232
4233 cs_tree_node_t *tn_p
4234 = cs_tree_get_node(cs_glob_tree,
4235 "thermophysical_models/porosities/porosity");
4236
4237 for (int z_id = 0; z_id < n_zones; z_id++) {
4238 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
4239
4240 if (z->type & CS_VOLUME_ZONE_POROSITY) {
4241 cs_tree_node_t *tn = _add_zone_id_test_attribute(tn_p, z->id);
4242 tn = cs_tree_get_node(tn, "model");
4243 const char *mdl = cs_tree_node_get_value_str(tn);
4244
4245 cs_glob_porous_model = CS_MAX(1, cs_glob_porous_model);
4246 if (mdl) {
4247 if (cs_gui_strcmp(mdl, "anisotropic"))
4248 cs_glob_porous_model = 2;
4249 else if (cs_gui_strcmp(mdl, "integral"))
4250 cs_glob_porous_model = 3;
4251 }
4252 }
4253 }
4254 }
4255
4256 /*-----------------------------------------------------------------------------
4257 * Get initial value from property markup.
4258 *
4259 * parameters:
4260 * property_name <-- name of the property
4261 * value --> new initial value of the property
4262 *----------------------------------------------------------------------------*/
4263
4264 void
cs_gui_properties_value(const char * property_name,double * value)4265 cs_gui_properties_value(const char *property_name,
4266 double *value)
4267 {
4268 cs_tree_node_t *tn = cs_tree_find_node(cs_glob_tree, "property");
4269 while (tn != NULL) {
4270 const char *name_tn = cs_tree_node_get_child_value_str(tn, "name");
4271 if (cs_gui_strcmp(name_tn, property_name))
4272 break;
4273 else
4274 tn = cs_tree_find_node_next(cs_glob_tree, tn, "property");
4275 }
4276 tn = cs_tree_get_node(tn, "initial_value");
4277
4278 cs_gui_node_get_real(tn, value);
4279 }
4280
4281 /*-----------------------------------------------------------------------------
4282 * Get value of property markup for fluid of given id
4283 *
4284 * parameters:
4285 * fluid_id <-- fluid index
4286 * property_name <-- name of the property
4287 * value --> new initial value of the property
4288 *----------------------------------------------------------------------------*/
4289
4290 void
cs_gui_properties_value_by_fluid_id(const int fluid_id,const char * property_name,double * value)4291 cs_gui_properties_value_by_fluid_id(const int fluid_id,
4292 const char *property_name,
4293 double *value)
4294 {
4295 cs_tree_node_t *tn = cs_tree_find_node(cs_glob_tree, "property");
4296 while (tn != NULL) {
4297 const char *name_tn = cs_tree_node_get_child_value_str(tn, "name");
4298 if (cs_gui_strcmp(name_tn, property_name))
4299 break;
4300 else
4301 tn = cs_tree_find_node_next(cs_glob_tree, tn, "property");
4302 }
4303
4304 char *label = NULL;
4305
4306 const char *prefix = "value_";
4307 size_t len = strlen(prefix) + 1;
4308 BFT_MALLOC(label, len+1, char);
4309
4310 sprintf(label, "%s%1i", prefix, fluid_id);
4311
4312 tn = cs_tree_get_node(tn, label);
4313 cs_gui_node_get_real(tn, value);
4314
4315 BFT_FREE(label);
4316 }
4317
4318 /*-----------------------------------------------------------------------------
4319 * Get value of reference fluid properties parameter.
4320 *
4321 * parameters:
4322 * name <-- parameter name
4323 * value --> parameter value
4324 *----------------------------------------------------------------------------*/
4325
4326 void
cs_gui_fluid_properties_value(const char * param,double * value)4327 cs_gui_fluid_properties_value(const char *param,
4328 double *value)
4329 {
4330 cs_tree_node_t *tn
4331 = cs_tree_get_node(cs_glob_tree, "physical_properties/fluid_properties");
4332 tn = cs_tree_get_node(tn, param);
4333
4334 cs_gui_node_get_real(tn, value);
4335 }
4336
4337 /*----------------------------------------------------------------------------
4338 * Get thermal scalar model.
4339 *
4340 * return:
4341 * value of itherm*10 + (temperature variant flag), or -1 if not defined
4342 *----------------------------------------------------------------------------*/
4343
4344 int
cs_gui_thermal_model(void)4345 cs_gui_thermal_model(void)
4346 {
4347 int test = -1;
4348
4349 const char *model = cs_gui_get_thermophysical_model("thermal_scalar");
4350
4351 if (model == NULL)
4352 return test;
4353
4354 if (cs_gui_strcmp(model, "off"))
4355 test = 0;
4356 else {
4357 if (cs_gui_strcmp(model, "enthalpy"))
4358 test = 20;
4359 else if (cs_gui_strcmp(model, "temperature_kelvin"))
4360 test = 11;
4361 else if (cs_gui_strcmp(model, "temperature_celsius"))
4362 test = 10;
4363 else if (cs_gui_strcmp(model, "potential_temperature"))
4364 test = 12;
4365 else if (cs_gui_strcmp(model, "liquid_potential_temperature"))
4366 test = 13;
4367 else if (cs_gui_strcmp(model, "total_energy"))
4368 test = 30;
4369 else
4370 bft_error(__FILE__, __LINE__, 0,
4371 _("Invalid thermal model: %s\n"), model);
4372 }
4373
4374 return test;
4375 }
4376
4377 /*----------------------------------------------------------------------------
4378 * Time moments definition
4379 *----------------------------------------------------------------------------*/
4380
4381 void
cs_gui_time_moments(void)4382 cs_gui_time_moments(void)
4383 {
4384 int imom = 1;
4385 int restart = cs_restart_present();
4386
4387 /* Loop on time average definitions */
4388
4389 const char path0[] = "/analysis_control/time_averages/time_average";
4390
4391 for (cs_tree_node_t *tn = cs_tree_get_node(cs_glob_tree, path0);
4392 tn != NULL;
4393 tn = cs_tree_node_get_next_of_name(tn), imom++) {
4394
4395 const char *restart_name;
4396 cs_time_moment_restart_t restart_mode = CS_TIME_MOMENT_RESTART_AUTO;
4397
4398 const int *v_i;
4399 const cs_real_t *v_r;
4400
4401 /* Older files used "label", now use "name", so try both */
4402
4403 const char *m_name = cs_tree_node_get_tag(tn, "name");
4404
4405 if (m_name == NULL) {
4406 m_name = cs_tree_node_get_tag(tn, "label");
4407 if (m_name == NULL) /* if neither found, force error case */
4408 m_name = cs_gui_node_get_tag(tn, "name");
4409 }
4410
4411 v_i = cs_tree_node_get_child_values_int(tn, "time_step_start");
4412 int nt_start = (v_i != NULL) ? v_i[0] : 0;
4413
4414 v_r = cs_tree_node_get_child_values_real(tn, "time_start");
4415 double t_start = (v_r != NULL) ? v_r[0] : -1;
4416
4417 /* test on restart */
4418
4419 if (restart != 0) {
4420 v_i = cs_tree_node_get_child_values_int(tn, "restart_from_time_average");
4421 int restart_id = (v_i != NULL) ? v_i[0] : -2;
4422 cs_time_moment_restart_options_by_id(restart_id,
4423 &restart_mode,
4424 &restart_name);
4425 }
4426
4427 int n_m_fields = cs_tree_get_node_count(tn, "var_prop");
4428
4429 int *m_f_id;
4430 BFT_MALLOC(m_f_id, n_m_fields*2, int);
4431 int *m_c_id = m_f_id + n_m_fields;
4432
4433 int j = 0;
4434 for (cs_tree_node_t *tn_vp = cs_tree_node_get_child(tn, "var_prop");
4435 tn_vp != NULL;
4436 tn_vp = cs_tree_node_get_next_of_name(tn_vp), j++) {
4437
4438 const char *f_name = cs_gui_node_get_tag(tn_vp, "name");
4439 v_i = cs_tree_node_get_child_values_int(tn_vp, "component");
4440 int idim = (v_i != NULL) ? v_i[0] : -1;
4441
4442 cs_field_t *f = cs_field_by_name_try(f_name);
4443
4444 if (f != NULL) {
4445 m_f_id[j] = f->id;
4446 m_c_id[j] = idim;
4447 }
4448
4449 else
4450 bft_error(__FILE__, __LINE__, 0,
4451 _("Time moment \"%s\"\n"
4452 "requires undefined field \"%s\"."),
4453 m_name, f_name);
4454
4455 }
4456
4457 cs_time_moment_define_by_field_ids(m_name,
4458 n_m_fields,
4459 m_f_id,
4460 m_c_id,
4461 CS_TIME_MOMENT_MEAN,
4462 nt_start,
4463 t_start,
4464 restart_mode,
4465 restart_name);
4466
4467 m_c_id = NULL;
4468 BFT_FREE(m_f_id);
4469
4470 }
4471
4472 #if _XML_DEBUG_
4473 bft_printf("==> %s\n", __func__);
4474 #endif
4475 }
4476
4477 /*-----------------------------------------------------------------------------
4478 * Set turbomachinery model
4479 *----------------------------------------------------------------------------*/
4480
4481 void
cs_gui_turbomachinery(void)4482 cs_gui_turbomachinery(void)
4483 {
4484 cs_turbomachinery_model_t model_type;
4485 bool coupled;
4486
4487 _turbomachinery_model(&model_type, &coupled);
4488
4489 cs_turbomachinery_set_model(model_type);
4490 }
4491
4492 /*-----------------------------------------------------------------------------
4493 * Set turbomachinery options.
4494 *----------------------------------------------------------------------------*/
4495
4496 void
cs_gui_turbomachinery_rotor(void)4497 cs_gui_turbomachinery_rotor(void)
4498 {
4499 cs_turbomachinery_model_t model_type;
4500 bool coupled;
4501
4502 _turbomachinery_model(&model_type, &coupled);
4503
4504 if (model_type != CS_TURBOMACHINERY_NONE) {
4505
4506 cs_tree_node_t *tn = NULL;
4507 cs_tree_node_t *tn2 = NULL;
4508 int n_rotors =
4509 cs_tree_get_node_count(cs_glob_tree,
4510 "/thermophysical_models/turbomachinery/rotor");
4511
4512 for (int rotor_id = 0; rotor_id < n_rotors; rotor_id++) {
4513
4514 double rotation_axis[3];
4515 double rotation_invariant[3];
4516 double rotation_velocity;
4517
4518 const char *cell_criteria;
4519
4520 rotation_axis[0] = _rotor_option(rotor_id, "axis_x");
4521 rotation_axis[1] = _rotor_option(rotor_id, "axis_y");
4522 rotation_axis[2] = _rotor_option(rotor_id, "axis_z");
4523
4524 rotation_invariant[0] = _rotor_option(rotor_id, "invariant_x");
4525 rotation_invariant[1] = _rotor_option(rotor_id, "invariant_y");
4526 rotation_invariant[2] = _rotor_option(rotor_id, "invariant_z");
4527
4528 tn = cs_tree_get_node(cs_glob_tree,
4529 "thermophysical_models/turbomachinery/rotor");
4530 int i;
4531 for (i = 1;
4532 tn != NULL && i < rotor_id + 1 ;
4533 i++) {
4534 tn = cs_tree_node_get_next_of_name(tn);
4535 }
4536 tn2 = tn;
4537 tn = cs_tree_get_node(tn, "velocity/value");
4538 cs_gui_node_get_real(tn, &rotation_velocity);
4539
4540 tn2 = cs_tree_get_node(tn2, "criteria");
4541 cell_criteria = cs_tree_node_get_value_str(tn2);
4542
4543 cs_turbomachinery_add_rotor(cell_criteria,
4544 rotation_velocity,
4545 rotation_axis,
4546 rotation_invariant);
4547
4548 }
4549
4550 int n_join = cs_tree_get_node_count(cs_glob_tree,
4551 "/thermophysical_models"
4552 "/turbomachinery/joining/face_joining");
4553
4554 for (int join_id = 0; join_id < n_join; join_id++) {
4555
4556 const char *selector_s = _get_rotor_face_joining("selector", join_id+1);
4557 const char *fraction_s = _get_rotor_face_joining("fraction", join_id+1);
4558 const char *plane_s = _get_rotor_face_joining("plane", join_id+1);
4559 const char *verbosity_s = _get_rotor_face_joining("verbosity", join_id+1);
4560 const char *visu_s = _get_rotor_face_joining("visualization", join_id+1);
4561
4562 double fraction = (fraction_s != NULL) ? atof(fraction_s) : 0.1;
4563 double plane = (plane_s != NULL) ? atof(plane_s) : 25.0;
4564 int verbosity = (verbosity_s != NULL) ? atoi(verbosity_s) : 0;
4565 int visualization = (visu_s != NULL) ? atoi(visu_s) : 0;
4566
4567 if (coupled == false)
4568 (void)cs_turbomachinery_join_add(selector_s,
4569 fraction,
4570 plane,
4571 verbosity,
4572 visualization);
4573 else
4574 (void)cs_turbomachinery_coupling_add(selector_s,
4575 fraction,
4576 verbosity);
4577
4578 }
4579
4580 }
4581 }
4582
4583 /*----------------------------------------------------------------------------
4584 * Logging output for MEI usage.
4585 *----------------------------------------------------------------------------*/
4586
4587 void
cs_gui_usage_log(void)4588 cs_gui_usage_log(void)
4589 {
4590 double mei_wtime = cs_gui_get_mei_times();
4591
4592 #if defined(HAVE_MPI)
4593
4594 if (cs_glob_n_ranks > 1) {
4595 double _wtime_loc = mei_wtime;
4596 MPI_Allreduce(&_wtime_loc, &mei_wtime, 1, MPI_DOUBLE, MPI_MAX,
4597 cs_glob_mpi_comm);
4598 }
4599
4600 #endif
4601
4602 if (mei_wtime > 0.0) {
4603 cs_log_printf(CS_LOG_PERFORMANCE,
4604 _("\nTime elapsed defining values using MEI: %12.5f\n"),
4605 mei_wtime);
4606 cs_log_printf(CS_LOG_PERFORMANCE, "\n");
4607 cs_log_separator(CS_LOG_PERFORMANCE);
4608 }
4609 }
4610
4611 /*----------------------------------------------------------------------------
4612 * Define user variables through the GUI.
4613 *----------------------------------------------------------------------------*/
4614
4615 void
cs_gui_user_variables(void)4616 cs_gui_user_variables(void)
4617 {
4618 int i = 0;
4619
4620 const char *t_scalar_name = NULL; /* thermal scalar name if present */
4621
4622 const char path_s[] = "additional_scalars/variable";
4623 cs_tree_node_t *tn_s = cs_tree_get_node(cs_glob_tree, path_s);
4624
4625 for (cs_tree_node_t *tn = tn_s;
4626 tn != NULL;
4627 tn = cs_tree_node_get_next_of_name(tn), i++) {
4628
4629 if (i == 0 && cs_glob_thermal_model->itherm != CS_THERMAL_MODEL_NONE) {
4630 const char path_t[] = "thermophysical_models/thermal_scalar/variable";
4631 t_scalar_name = cs_tree_node_get_tag
4632 (cs_tree_get_node(cs_glob_tree, path_t), "name");
4633 }
4634
4635 const char *name = cs_gui_node_get_tag(tn, "name");
4636
4637 const char *variance_name = cs_tree_node_get_child_value_str(tn, "variance");
4638
4639 /* In case of variance, check for presence of matching field
4640 in thermal and user scalars */
4641
4642 if (variance_name != NULL) {
4643
4644 bool found = false;
4645 if (t_scalar_name != NULL) {
4646 if (strcmp(t_scalar_name, variance_name) == 0)
4647 found = true;
4648 }
4649 for (cs_tree_node_t *tn_c = tn_s;
4650 tn_c != NULL && found == false;
4651 tn_c = cs_tree_node_get_next_of_name(tn_c), i++) {
4652 const char *cmp_name = cs_tree_node_get_tag(tn_c, "name");
4653 if (cmp_name != NULL) {
4654 if (strcmp(cmp_name, variance_name) == 0)
4655 found = true;
4656 }
4657 }
4658
4659 if (found)
4660 cs_parameters_add_variable_variance(name, variance_name);
4661
4662 }
4663
4664 /* If not a variance, we have a regular variable */
4665
4666 else
4667 cs_parameters_add_variable(name, 1);
4668 }
4669 }
4670
4671 /*----------------------------------------------------------------------------
4672 * Define user arrays through the GUI.
4673 *----------------------------------------------------------------------------*/
4674
4675 void
cs_gui_user_arrays(void)4676 cs_gui_user_arrays(void)
4677 {
4678 const char path_s[] = "additional_scalars/users/property";
4679 cs_tree_node_t *tn_s = cs_tree_get_node(cs_glob_tree, path_s);
4680
4681 for (cs_tree_node_t *tn = tn_s;
4682 tn != NULL;
4683 tn = cs_tree_node_get_next_of_name(tn)) {
4684
4685 const char *name = cs_gui_node_get_tag(tn, "name");
4686
4687 int array_dim = 1;
4688 cs_tree_node_t *dtn = cs_tree_get_node(tn, "dimension");
4689 cs_gui_node_get_int(dtn, &array_dim);
4690
4691 const char *location_name = cs_gui_node_get_tag(tn, "support");
4692
4693 if (strcmp(location_name, "cells") == 0)
4694 cs_parameters_add_property(name, array_dim, CS_MESH_LOCATION_CELLS);
4695
4696 else if (strcmp(location_name, "internal") == 0)
4697 cs_parameters_add_property(name,
4698 array_dim,
4699 CS_MESH_LOCATION_INTERIOR_FACES);
4700
4701 else if (strcmp(location_name, "boundary") == 0)
4702 cs_parameters_add_property(name,
4703 array_dim,
4704 CS_MESH_LOCATION_BOUNDARY_FACES);
4705
4706 else if (strcmp(location_name, "vertices") == 0)
4707 cs_parameters_add_property(name,
4708 array_dim,
4709 CS_MESH_LOCATION_VERTICES);
4710
4711 }
4712 }
4713
4714 /*----------------------------------------------------------------------------
4715 * Define volume and boundary zones through the GUI.
4716 *----------------------------------------------------------------------------*/
4717
4718 void
cs_gui_zones(void)4719 cs_gui_zones(void)
4720 {
4721 /* Ensure zones ordering for safety (should be removed in the future)*/
4722
4723 _ensure_zones_order();
4724
4725 int id = 0;
4726
4727 const char default_criteria[] = "all[]";
4728
4729 /* Volume zones */
4730 /*------------- */
4731
4732 cs_tree_node_t *tn_vc = cs_tree_get_node(cs_glob_tree,
4733 "solution_domain/volumic_conditions");
4734
4735 const int n_v_zones = cs_tree_get_node_count(tn_vc, "zone");
4736
4737 /* Build ordering array to check zones are defined in increasing id order */
4738
4739 cs_lnum_t *order = NULL, *z_ids = NULL;
4740 BFT_MALLOC(order, n_v_zones, cs_lnum_t);
4741 BFT_MALLOC(z_ids, n_v_zones, cs_lnum_t);
4742
4743 /* Loop on volume condition zones */
4744
4745 id = 0;
4746 for (cs_tree_node_t *tn = cs_tree_node_get_child(tn_vc, "zone");
4747 tn != NULL;
4748 tn = cs_tree_node_get_next_of_name(tn), id++) {
4749 z_ids[id] = _v_zone_t_id(tn, id);
4750 }
4751
4752 assert(id == n_v_zones);
4753
4754 cs_order_lnum_allocated(NULL, z_ids, order, n_v_zones);
4755
4756 /* Now loop on zones in id order */
4757
4758 for (int i = 0; i < n_v_zones; i++) {
4759
4760 int type_flag = 0, z_id = z_ids[order[i]];
4761
4762 cs_tree_node_t *tn = _v_zone_node_by_id(tn_vc, z_id);
4763
4764 /* zone name */
4765
4766 const char *name = cs_tree_node_get_tag(tn, "label");
4767
4768 /* location criteria */
4769
4770 const char *_criteria = cs_tree_node_get_value_str(tn);
4771 const char *criteria = (_criteria != NULL) ? _criteria : default_criteria;
4772
4773 /* Check for initialization */
4774
4775 if (_zone_is_type(tn, "initialization"))
4776 type_flag = type_flag | CS_VOLUME_ZONE_INITIALIZATION;
4777
4778 /* Check for porosity */
4779
4780 if (_zone_is_type(tn, "porosity"))
4781 type_flag = type_flag | CS_VOLUME_ZONE_POROSITY;
4782
4783 /* Check for head losses */
4784
4785 if (_zone_is_type(tn, "head_losses"))
4786 type_flag = type_flag | CS_VOLUME_ZONE_HEAD_LOSS;
4787
4788 /* Check for source terms */
4789
4790 if (_zone_is_type(tn, "momentum_source_term"))
4791 type_flag = type_flag | CS_VOLUME_ZONE_SOURCE_TERM;
4792 if (_zone_is_type(tn, "scalar_source_term"))
4793 type_flag = type_flag | CS_VOLUME_ZONE_SOURCE_TERM;
4794 if (_zone_is_type(tn, "thermal_source_term"))
4795 type_flag = type_flag | CS_VOLUME_ZONE_SOURCE_TERM;
4796
4797 /* Check for solid zones */
4798
4799 if (_zone_is_type(tn, "solid"))
4800 type_flag = type_flag | CS_VOLUME_ZONE_SOLID;
4801
4802 /* Check if zone is used to define variable physical properties */
4803 if (_zone_is_type(tn, "physical_properties"))
4804 type_flag = type_flag | CS_VOLUME_ZONE_PHYSICAL_PROPERTIES;
4805 else {
4806 /* FIXME: zone all_cells is used for physical properties no matter what.
4807 * Will be changed in v7.1 for a cleaner definition
4808 */
4809 if (cs_gui_strcmp(name, "all_cells"))
4810 type_flag = type_flag | CS_VOLUME_ZONE_PHYSICAL_PROPERTIES;
4811 }
4812
4813 /* Finally, define zone */
4814
4815 cs_volume_zone_define(name, criteria, type_flag);
4816 }
4817
4818 BFT_FREE(order);
4819 BFT_FREE(z_ids);
4820
4821 /* Boundary zones */
4822 /*--------------- */
4823
4824 /* Loop on boundary condition zones */
4825
4826 cs_tree_node_t *tn_bc = cs_tree_get_node(cs_glob_tree,
4827 "boundary_conditions");
4828
4829 id = 0;
4830 for (cs_tree_node_t *tn = cs_tree_node_get_child(tn_bc, "boundary");
4831 tn != NULL;
4832 tn = cs_tree_node_get_next_of_name(tn), id++) {
4833
4834 /* Zone id in tree; note that the name tag for boundary zones actually
4835 defines an integer (1 to n). This tag should be removed in the
4836 future to only use the zone label (the actual zone name). */
4837
4838 const char *id_s = cs_tree_node_get_tag(tn, "name");
4839 if (id_s != NULL) {
4840 int z_t_id = atoi(id_s);
4841 if (z_t_id != id + 1)
4842 bft_printf(_("\n"
4843 " Warning: noncontiguous %s zone ids in XML:\n"
4844 " zone with index %d has id %d.\n"),
4845 tn->name, id, z_t_id);
4846 }
4847
4848 /* Zone name */
4849
4850 const char *name = cs_tree_node_get_tag(tn, "label");
4851
4852 /* location criteria */
4853
4854 const char *_criteria = cs_tree_node_get_value_str(tn);
4855 const char *criteria = (_criteria != NULL) ? _criteria : default_criteria;
4856
4857 int type_flag = 0;
4858
4859 /* Define zone */
4860
4861 cs_boundary_zone_define(name, criteria, type_flag);
4862 }
4863
4864 }
4865
4866 /*----------------------------------------------------------------------------
4867 * Define balance by zone through the GUI.
4868 *----------------------------------------------------------------------------*/
4869
4870 void
cs_gui_balance_by_zone(void)4871 cs_gui_balance_by_zone(void)
4872 {
4873 const char path0[] = "/analysis_control/scalar_balances/scalar_balance";
4874
4875 for (cs_tree_node_t *tn = cs_tree_get_node(cs_glob_tree, path0);
4876 tn != NULL;
4877 tn = cs_tree_node_get_next_of_name(tn)) {
4878
4879 const char _default_criteria[] = "all[]";
4880
4881 const char *criteria = cs_tree_node_get_child_value_str(tn, "criteria");
4882 if (criteria == NULL) criteria = _default_criteria;
4883
4884 for (cs_tree_node_t *tn_v = cs_tree_node_get_child(tn, "var_prop");
4885 tn_v != NULL;
4886 tn_v = cs_tree_node_get_next_of_name(tn_v)) {
4887
4888 const char *name = cs_gui_node_get_tag(tn_v, "name");
4889 cs_balance_by_zone(criteria, name);
4890
4891 }
4892 }
4893 }
4894
4895 /*----------------------------------------------------------------------------
4896 * Define pressure drop through the GUI.
4897 *----------------------------------------------------------------------------*/
4898
4899 void
cs_gui_pressure_drop_by_zone(void)4900 cs_gui_pressure_drop_by_zone(void)
4901 {
4902 const char path0[] = "/analysis_control/scalar_balances/pressure_drop";
4903
4904 for (cs_tree_node_t *tn = cs_tree_get_node(cs_glob_tree, path0);
4905 tn != NULL;
4906 tn = cs_tree_node_get_next_of_name(tn)) {
4907
4908 const char _default_criteria[] = "all[]";
4909
4910 const char *criteria = cs_tree_node_get_child_value_str(tn, "criteria");
4911 if (criteria == NULL) criteria = _default_criteria;
4912
4913 cs_pressure_drop_by_zone(criteria);
4914 }
4915 }
4916
4917 /*----------------------------------------------------------------------------
4918 * Define fans through the GUI.
4919 *----------------------------------------------------------------------------*/
4920
4921 void
cs_gui_define_fans(void)4922 cs_gui_define_fans(void)
4923 {
4924 const char path0[] = "thermophysical_models/fans/fan";
4925
4926 for (cs_tree_node_t *tn = cs_tree_get_node(cs_glob_tree, path0);
4927 tn != NULL;
4928 tn = cs_tree_node_get_next_of_name(tn)) {
4929
4930 const int *v_i;
4931 const cs_real_t *v_r;
4932
4933 const char *i_axis_s[] = {"inlet_axis_x", "inlet_axis_y", "inlet_axis_z"};
4934 const char *o_axis_s[] = {"outlet_axis_x", "outlet_axis_y", "outlet_axis_z"};
4935 const char *p_coeff_s[]
4936 = {"curve_coeffs_x", "curve_coeffs_y", "curve_coeffs_z"};
4937
4938 v_i = cs_tree_node_get_child_values_int(tn, "mesh_dimension");
4939 int dim = (v_i != NULL) ? v_i[0] : 3;
4940
4941 cs_real_t inlet_axis_coords[3] = {0, 0, 0};
4942 cs_real_t outlet_axis_coords[3] = {0.1, 0, 0};
4943 cs_real_t pressure_curve_coeffs[3] = {0.6, -0.1, -0.05};
4944
4945 for (int i = 0; i < 3; i++) {
4946 v_r = cs_tree_node_get_child_values_real(tn, i_axis_s[i]);
4947 if (v_r != NULL) inlet_axis_coords[i] = v_r[0];
4948 }
4949 for (int i = 0; i < 3; i++) {
4950 v_r = cs_tree_node_get_child_values_real(tn, o_axis_s[i]);
4951 if (v_r != NULL) outlet_axis_coords[i] = v_r[0];
4952 }
4953
4954 v_r = cs_tree_node_get_child_values_real(tn, "fan_radius");
4955 cs_real_t fan_radius = (v_r != NULL) ? v_r[0] : 0.7;
4956
4957 v_r = cs_tree_node_get_child_values_real(tn, "blades_radius");
4958 cs_real_t blades_radius = (v_r != NULL) ? v_r[0] : 0.5;
4959
4960 v_r = cs_tree_node_get_child_values_real(tn, "hub_radius");
4961 cs_real_t hub_radius = (v_r != NULL) ? v_r[0] : 0.1;
4962
4963 v_r = cs_tree_node_get_child_values_real(tn, "axial_torque");
4964 cs_real_t axial_torque = (v_r != NULL) ? v_r[0] : 0.01;
4965
4966 for (int i = 0; i < 3; i++) {
4967 v_r = cs_tree_node_get_child_values_real(tn, p_coeff_s[i]);
4968 if (v_r != NULL) pressure_curve_coeffs[i] = v_r[0];
4969 }
4970
4971 cs_fan_define(dim, /* fan (mesh) dimension (2D or 3D) */
4972 0, /* mode */
4973 inlet_axis_coords,
4974 outlet_axis_coords,
4975 fan_radius,
4976 blades_radius,
4977 hub_radius,
4978 pressure_curve_coeffs,
4979 axial_torque);
4980 }
4981 }
4982
4983 /*----------------------------------------------------------------------------
4984 * Define error estimator through the GUI.
4985 *----------------------------------------------------------------------------*/
4986
4987 void
cs_gui_error_estimator(int * iescal,int * iespre,int * iesder,int * iescor,int * iestot)4988 cs_gui_error_estimator(int *iescal,
4989 int *iespre,
4990 int *iesder,
4991 int *iescor,
4992 int *iestot)
4993 {
4994 cs_tree_node_t *tn_ee
4995 = cs_tree_get_node(cs_glob_tree, "analysis_control/error_estimator");
4996
4997 cs_tree_node_t *tn = cs_tree_get_node(tn_ee, "Correction/model");
4998
4999 const char *result = cs_tree_node_get_value_str(tn);
5000
5001 if (cs_gui_strcmp(result, "1"))
5002 iescal[*iescor -1] = 1;
5003 else if (cs_gui_strcmp(result, "2"))
5004 iescal[*iescor -1] = 2;
5005 else
5006 iescal[*iescor -1] = 0;
5007
5008 tn = cs_tree_get_node(tn_ee, "Drift/model");
5009
5010 result = cs_tree_node_get_value_str(tn);
5011
5012 if (cs_gui_strcmp(result, "1"))
5013 iescal[*iesder -1] = 1;
5014 else if (cs_gui_strcmp(result, "2"))
5015 iescal[*iesder -1] = 2;
5016 else
5017 iescal[*iesder -1] = 0;
5018
5019 tn = cs_tree_get_node(tn_ee, "Prediction/model");
5020
5021 result = cs_tree_node_get_value_str(tn);
5022
5023 if (cs_gui_strcmp(result, "1"))
5024 iescal[*iespre -1] = 1;
5025 else if (cs_gui_strcmp(result, "2"))
5026 iescal[*iespre -1] = 2;
5027 else
5028 iescal[*iespre -1] = 0;
5029
5030 tn = cs_tree_get_node(tn_ee, "Total/model");
5031
5032 result = cs_tree_node_get_value_str(tn);
5033
5034 if (cs_gui_strcmp(result, "1"))
5035 iescal[*iestot -1] = 1;
5036 else if (cs_gui_strcmp(result, "2"))
5037 iescal[*iestot -1] = 2;
5038 else
5039 iescal[*iestot -1] = 0;
5040 }
5041
5042 /*----------------------------------------------------------------------------
5043 * Define internal coupling through the GUI.
5044 *----------------------------------------------------------------------------*/
5045
5046 void
cs_gui_internal_coupling(void)5047 cs_gui_internal_coupling(void)
5048 {
5049 int n_zones = cs_volume_zone_n_zones();
5050
5051 int n_solid_zones = 0;
5052 for (int i = 0; i < n_zones; i++) {
5053 const cs_zone_t *z = cs_volume_zone_by_id(i);
5054 if (z->type & CS_VOLUME_ZONE_SOLID)
5055 n_solid_zones += 1;
5056 }
5057
5058 if (n_solid_zones < 1)
5059 return;
5060
5061 cs_tree_node_t *node_int_cpl
5062 = cs_tree_get_node(cs_glob_tree, "thermophysical_models/internal_coupling");
5063
5064 if (node_int_cpl != NULL) {
5065
5066 int j = 0;
5067 int *z_ids;
5068 BFT_MALLOC(z_ids, n_solid_zones, int);
5069
5070 for (int i = 0; i < n_zones; i++) {
5071 const cs_zone_t *z = cs_volume_zone_by_id(i);
5072 if (z->type & CS_VOLUME_ZONE_SOLID)
5073 z_ids[j++] = z->id;
5074 }
5075
5076 int coupling_id = cs_internal_coupling_n_couplings();
5077
5078 cs_internal_coupling_add_volume_zones(n_solid_zones, z_ids);
5079 BFT_FREE(z_ids);
5080
5081 {
5082 cs_internal_coupling_t *cpl = cs_internal_coupling_by_id(coupling_id);
5083
5084 char i_name[64], e_name[64];
5085 snprintf(i_name, 63, "auto:internal_coupling_%d_fluid", cpl->id);
5086 i_name[63] = '\0';
5087 snprintf(e_name, 63, "auto:internal_coupling_%d_solid", cpl->id);
5088 e_name[63] = '\0';
5089
5090 cs_internal_coupling_add_boundary_groups(cpl, i_name, e_name);
5091 }
5092
5093 if (n_solid_zones > 0) {
5094 cs_tree_node_t *ns
5095 = cs_tree_node_get_child(node_int_cpl, "coupled_scalars");
5096 /* Add the coupled scalars defined in the GUI */
5097 for (cs_tree_node_t *tn = cs_tree_node_get_child(ns, "scalar");
5098 tn != NULL;
5099 tn = cs_tree_node_get_next_of_name(tn)) {
5100
5101 const char *scalar_name = cs_tree_node_get_tag(tn, "name");
5102 int f_id = cs_field_id_by_name(scalar_name);
5103 if (f_id < 0)
5104 bft_error(__FILE__, __LINE__, 0,
5105 _("Scalar %s does not exist!.\n"), scalar_name);
5106
5107 cs_internal_coupling_add_entity(f_id);
5108 }
5109 }
5110 }
5111 }
5112
5113 /*----------------------------------------------------------------------------*/
5114
5115 END_C_DECLS
5116