1 /*============================================================================
2 * Physical properties computation and management.
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 <assert.h>
34 #include <math.h>
35 #include <stdarg.h>
36 #include <stdlib.h>
37 #include <stdio.h>
38 #include <string.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *----------------------------------------------------------------------------*/
43
44 #include "bft_error.h"
45 #include "bft_mem.h"
46 #include "bft_printf.h"
47
48 #include "cs_property.h"
49
50 /*----------------------------------------------------------------------------
51 * Header for the current file
52 *----------------------------------------------------------------------------*/
53
54 #include "cs_physical_properties.h"
55 #if defined(HAVE_EOS)
56 #include "cs_eos.hxx"
57 #endif
58
59 #if defined(HAVE_FREESTEAM)
60 #include <freesteam/steam_ph.h>
61 #include <freesteam/steam_pT.h>
62 #include <freesteam/steam_ps.h>
63 #include <freesteam/steam_pu.h>
64 #include <freesteam/steam_pv.h>
65 #include <freesteam/steam_Ts.h>
66 #include <freesteam/steam_Tx.h>
67
68 #include <freesteam/region1.h>
69 #include <freesteam/region2.h>
70 #include <freesteam/region3.h>
71 #include <freesteam/region4.h>
72 #endif
73
74 #if defined(HAVE_COOLPROP)
75 #include "cs_coolprop.hxx"
76 #endif
77
78 /*----------------------------------------------------------------------------*/
79
80 BEGIN_C_DECLS
81
82 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
83
84 /*============================================================================
85 * Local macro definitions
86 *============================================================================*/
87
88 /* Directory name separator
89 (historically, '/' for Unix/Linux, '\' for Windows, ':' for Mac
90 but '/' should work for all on modern systems) */
91
92 #define DIR_SEPARATOR '/'
93
94 /*============================================================================
95 * Type definitions
96 *============================================================================*/
97
98 /* Thermal table structure */
99
100 typedef struct {
101
102 char *material; /* material choice (water, ...) */
103 char *method; /* method choice
104 (cathare, thetis, freesteam, ...) */
105 int type; /* 0 for user
106 * 1 for freesteam
107 * 2 for eos
108 * 3 for coolprop
109 */
110
111 cs_phys_prop_thermo_plane_type_t thermo_plane;
112
113 int temp_scale; /* temperature scale if needed
114 * 1 for kelvin
115 * 2 for Celsius */
116
117 } cs_thermal_table_t;
118
119 /*----------------------------------------------------------------------------
120 * Function pointer types
121 *----------------------------------------------------------------------------*/
122
123 typedef void
124 (cs_eos_create_t)(char *EOSMethod,
125 char *EOSRef);
126
127 typedef void
128 (cs_eos_destroy_t)(void);
129
130 typedef void
131 (cs_phys_prop_eos_t)(cs_phys_prop_thermo_plane_type_t thermo_plane,
132 cs_phys_prop_type_t property,
133 const cs_lnum_t n_vals,
134 double var1[],
135 double var2[],
136 cs_real_t val[]);
137
138 typedef void
139 (cs_phys_prop_coolprop_t)(const char *coolprop_material,
140 const char *coolprop_backend,
141 cs_phys_prop_thermo_plane_type_t thermo_plane,
142 cs_phys_prop_type_t property,
143 const cs_lnum_t n_vals,
144 const cs_real_t var1[],
145 const cs_real_t var2[],
146 cs_real_t val[]);
147
148 /*============================================================================
149 * Static global variables
150 *============================================================================*/
151
152 cs_thermal_table_t *cs_glob_thermal_table = NULL;
153
154 #if defined(HAVE_DLOPEN) && defined(HAVE_EOS)
155
156 static void *_cs_eos_dl_lib = NULL;
157 static cs_eos_create_t *_cs_eos_create = NULL;
158 static cs_eos_destroy_t *_cs_eos_destroy = NULL;
159 static cs_phys_prop_eos_t *_cs_phys_prop_eos = NULL;
160
161 #endif
162
163 #if defined(HAVE_COOLPROP)
164
165 static cs_phys_prop_coolprop_t *_cs_phys_prop_coolprop = NULL;
166 #if defined(HAVE_PLUGINS)
167 static void *_cs_coolprop_dl_lib = NULL;
168 #endif
169
170 #endif
171
172
173 /*============================================================================
174 * Private function definitions
175 *============================================================================*/
176
177 /*----------------------------------------------------------------------------
178 * Create an empty thermal_table structure
179 *----------------------------------------------------------------------------*/
180
181 static cs_thermal_table_t *
_thermal_table_create(void)182 _thermal_table_create(void)
183 {
184 cs_thermal_table_t *tt = NULL;
185
186 BFT_MALLOC(tt, 1, cs_thermal_table_t);
187
188 tt->material = NULL;
189 tt->method = NULL;
190 tt->type = 0;
191 tt->temp_scale = 0;
192 tt->thermo_plane = CS_PHYS_PROP_PLANE_PH;
193
194 return tt;
195 }
196
197 /*----------------------------------------------------------------------------*/
198 /*!
199 * \brief Get xdef of a property on a given zone.
200 *
201 * \param[in] pty pointer to cs_property_t
202 * \param[in] zname name of the zone. Can be NULL for 'all_cells'
203 *
204 * \return pointer to corresponding cs_xdef_t
205 *
206 */
207 /*----------------------------------------------------------------------------*/
208
209 static cs_xdef_t *
_get_property_def_on_zone(const cs_property_t * pty,const char * zname)210 _get_property_def_on_zone(const cs_property_t *pty,
211 const char *zname)
212 {
213 cs_xdef_t *def = NULL;
214
215 const int z_id = cs_get_vol_zone_id(zname);
216 for (int i = 0; i < pty->n_definitions; i++) {
217 if (pty->defs[i]->z_id == z_id) {
218 def = pty->defs[i];
219 break;
220 }
221 }
222
223 return def;
224 }
225
226 /*----------------------------------------------------------------------------*/
227 /*!
228 * \brief update values of a definition
229 *
230 * \param[in] def pointer to cs_xdef_t structure
231 * \param[in] new_vals array of new values
232 *
233 */
234 /*----------------------------------------------------------------------------*/
235
236 static void
_update_def_values(cs_xdef_t * def,const cs_real_t * new_vals)237 _update_def_values(cs_xdef_t *def,
238 const cs_real_t *new_vals)
239 {
240 cs_real_t *_context = (cs_real_t *)def->context;
241
242 for (int i = 0; i < def->dim; i++)
243 _context[i] = new_vals[i];
244 }
245
246 /*----------------------------------------------------------------------------*/
247 /*!
248 * \brief Create a new physical property using the cs_property_t structure.
249 *
250 * \param[in] name name of the property
251 * \param[in] dim dimension of the property (1->scalar, 3->vector,..)
252 * \param[in] refval Reference value
253 *
254 * \return pointer to the newly created cs_property_t structure.
255 *
256 */
257 /*----------------------------------------------------------------------------*/
258
259 static cs_property_t *
_physical_property_create(const char * name,const int dim,const cs_real_t refval)260 _physical_property_create(const char *name,
261 const int dim,
262 const cs_real_t refval)
263 {
264 cs_property_t *pty = cs_property_by_name(name);
265 if (pty != NULL)
266 bft_error(__FILE__, __LINE__, 0,
267 _("Error: property '%s' is already defined.\n"),
268 name);
269
270 if (dim == 1)
271 pty = cs_property_add(name, CS_PROPERTY_ISO);
272 else if (dim == 3)
273 pty = cs_property_add(name, CS_PROPERTY_ORTHO);
274 else if (dim == 6)
275 pty = cs_property_add(name, CS_PROPERTY_ANISO_SYM);
276 else if (dim == 9)
277 pty = cs_property_add(name, CS_PROPERTY_ANISO);
278 else
279 bft_error(__FILE__, __LINE__, 0,
280 _("Error: for property '%s', dimension %d not supported.\n"),
281 name, dim);
282
283 cs_property_set_reference_value(pty, refval);
284
285 return pty;
286 }
287
288 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
289
290 /*=============================================================================
291 * Public function definitions
292 *============================================================================*/
293
294 /*----------------------------------------------------------------------------*/
295 /*!
296 * \brief Define thermal table.
297 */
298 /*----------------------------------------------------------------------------*/
299
300 void
cs_thermal_table_set(const char * material,const char * method,const char * reference,cs_phys_prop_thermo_plane_type_t thermo_plane,int temp_scale)301 cs_thermal_table_set(const char *material,
302 const char *method,
303 const char *reference,
304 cs_phys_prop_thermo_plane_type_t thermo_plane,
305 int temp_scale)
306 {
307 if (cs_glob_thermal_table == NULL)
308 cs_glob_thermal_table = _thermal_table_create();
309
310 BFT_MALLOC(cs_glob_thermal_table->material, strlen(material) +1, char);
311 strcpy(cs_glob_thermal_table->material, material);
312
313 if (strcmp(method, "freesteam") == 0 ||
314 strcmp(material, "user_material") == 0) {
315 BFT_MALLOC(cs_glob_thermal_table->method, strlen(method) +1, char);
316 if (strcmp(method, "freesteam") == 0)
317 cs_glob_thermal_table->type = 1;
318 else
319 cs_glob_thermal_table->type = 0;
320 }
321 else if (strcmp(method, "CoolProp") == 0) {
322 BFT_MALLOC(cs_glob_thermal_table->method, strlen(method) +1, char);
323 cs_glob_thermal_table->type = 3;
324 #if defined(HAVE_COOLPROP)
325 #if defined(HAVE_PLUGINS)
326 {
327 /* Open from shared library */
328 _cs_coolprop_dl_lib = cs_base_dlopen_plugin("cs_coolprop");
329
330 /* Load symbols from shared library */
331
332 /* Function pointers need to be double-cast so as to first convert
333 a (void *) type to a memory address and then convert it back to the
334 original type. Otherwise, the compiler may issue a warning.
335 This is a valid ISO C construction. */
336
337 _cs_phys_prop_coolprop = (cs_phys_prop_coolprop_t *) (intptr_t)
338 cs_base_get_dl_function_pointer(_cs_coolprop_dl_lib,
339 "cs_phys_prop_coolprop",
340 true);
341 }
342 #else
343 _cs_phys_prop_coolprop = (cs_phys_prop_coolprop_t *) (intptr_t)
344 cs_phys_prop_coolprop;
345 #endif
346 #endif
347 }
348 else {
349 BFT_MALLOC(cs_glob_thermal_table->method, strlen(method) +5, char);
350 strcpy(cs_glob_thermal_table->method, "EOS_");
351 strcat(cs_glob_thermal_table->method, method);
352 cs_glob_thermal_table->type = 2;
353 #if defined(HAVE_EOS)
354 {
355 const char _reference_default[] = "";
356 const char *_reference = (reference != NULL) ? reference : _reference_default;
357
358 /* Open from shared library */
359 _cs_eos_dl_lib = cs_base_dlopen_plugin("cs_eos");
360
361 /* Function pointers need to be double-casted so as to first convert
362 a (void *) type to a memory address and then convert it back to the
363 original type. Otherwise, the compiler may issue a warning.
364 This is a valid ISO C construction. */
365
366 _cs_eos_create = (cs_eos_create_t *) (intptr_t)
367 cs_base_get_dl_function_pointer(_cs_eos_dl_lib, "cs_eos_create", true);
368 _cs_eos_destroy = (cs_eos_destroy_t *) (intptr_t)
369 cs_base_get_dl_function_pointer(_cs_eos_dl_lib, "cs_eos_destroy", true);
370 _cs_phys_prop_eos = (cs_phys_prop_eos_t *) (intptr_t)
371 cs_base_get_dl_function_pointer(_cs_eos_dl_lib, "cs_phys_prop_eos", true);
372
373 _cs_eos_create(cs_glob_thermal_table->method,
374 _reference);
375 }
376 #else
377 CS_UNUSED(reference);
378 #endif
379 }
380 cs_glob_thermal_table->thermo_plane = thermo_plane;
381 cs_glob_thermal_table->temp_scale = temp_scale;
382 }
383
384 /*----------------------------------------------------------------------------*/
385 /*!
386 * \brief finalize thermal table.
387 */
388 /*----------------------------------------------------------------------------*/
389
390 void
cs_thermal_table_finalize(void)391 cs_thermal_table_finalize(void)
392 {
393 if (cs_glob_thermal_table != NULL) {
394 #if defined(HAVE_EOS)
395 if (cs_glob_thermal_table->type == 2) {
396 _cs_eos_destroy();
397 cs_base_dlclose("cs_eos", _cs_eos_dl_lib);
398 _cs_eos_create = NULL;
399 _cs_eos_destroy = NULL;
400 _cs_phys_prop_eos = NULL;
401 }
402 #endif
403 #if defined(HAVE_COOLPROP) && defined(HAVE_PLUGINS)
404 if (cs_glob_thermal_table->type == 3) {
405 cs_base_dlclose("cs_coolprop", _cs_coolprop_dl_lib);
406 _cs_phys_prop_coolprop = NULL;
407 }
408 #endif
409 BFT_FREE(cs_glob_thermal_table->material);
410 BFT_FREE(cs_glob_thermal_table->method);
411 BFT_FREE(cs_glob_thermal_table);
412 }
413 }
414
415 /*----------------------------------------------------------------------------*/
416 /*!
417 * \brief Compute a physical property.
418 *
419 * For values var1 and var2, we can use a stride so that accesses for a given
420 * element with id i will be of the form: var[i*stride]; this allows regular
421 * access with stride 1, and access to constant variables stored as a
422 * single-valued array with a stride of 0.
423 *
424 * \param[in] property property queried
425 * \param[in] n_vals number of values
426 * \param[in] var1_stride stride between successive values of var1
427 * \param[in] var2_stride stride between successive values of var2
428 * \param[in] var1 values on first plane axis
429 * \param[in] var2 values on second plane axis
430 * \param[out] val resulting property values
431 */
432 /*----------------------------------------------------------------------------*/
433
434 void
cs_phys_prop_compute(cs_phys_prop_type_t property,cs_lnum_t n_vals,cs_lnum_t var1_stride,cs_lnum_t var2_stride,const cs_real_t var1[],const cs_real_t var2[],cs_real_t val[])435 cs_phys_prop_compute(cs_phys_prop_type_t property,
436 cs_lnum_t n_vals,
437 cs_lnum_t var1_stride,
438 cs_lnum_t var2_stride,
439 const cs_real_t var1[],
440 const cs_real_t var2[],
441 cs_real_t val[])
442 {
443 cs_lnum_t _n_vals = n_vals;
444 cs_real_t _var2_c_single[1];
445 cs_real_t *_var1_c = NULL, *_var2_c = NULL;
446 const cs_real_t *var1_c = var1, *var2_c = var2;
447
448 if (n_vals < 1)
449 return;
450
451 /* Adapt to different strides to optimize for constant arrays */
452
453 if (var1_stride == 0 && var2_stride == 0)
454 _n_vals = 1;
455
456 if (var1_stride == 0 && n_vals > 1) {
457 BFT_MALLOC(_var1_c, n_vals, cs_real_t);
458 for (cs_lnum_t ii = 0; ii < n_vals; ii++)
459 _var1_c[ii] = var1[0];
460 var1_c = _var1_c;
461 }
462
463 if (cs_glob_thermal_table->temp_scale == 2) {
464 if (_n_vals == 1) {
465 _var2_c_single[0] = var2[0] + 273.15;
466 var2_c = _var2_c_single;
467 }
468 else {
469 BFT_MALLOC(_var2_c, n_vals, cs_real_t);
470 for (cs_lnum_t ii = 0; ii < n_vals; ii++)
471 _var2_c[ii] = var2[ii*var2_stride] + 273.15;
472 var2_c = _var2_c;
473 }
474 }
475 else {
476 if (var2_stride == 0 && n_vals > 1) {
477 BFT_MALLOC(_var2_c, n_vals, cs_real_t);
478 for (cs_lnum_t ii = 0; ii < n_vals; ii++)
479 _var2_c[ii] = var2[0];
480 var2_c = _var2_c;
481 }
482 }
483
484 /* Compute proper */
485
486 if (cs_glob_thermal_table->type == 1) {
487 cs_phys_prop_freesteam(cs_glob_thermal_table->thermo_plane,
488 property,
489 _n_vals,
490 var1_c,
491 var2_c,
492 val);
493 }
494 #if defined(HAVE_EOS)
495 else if (cs_glob_thermal_table->type == 2) {
496 _cs_phys_prop_eos(cs_glob_thermal_table->thermo_plane,
497 property,
498 _n_vals,
499 var1_c,
500 var2_c,
501 val);
502 }
503 #endif
504 #if defined(HAVE_COOLPROP)
505 else if (cs_glob_thermal_table->type == 3) {
506 _cs_phys_prop_coolprop(cs_glob_thermal_table->material,
507 "HEOS",
508 cs_glob_thermal_table->thermo_plane,
509 property,
510 _n_vals,
511 var1_c,
512 var2_c,
513 val);
514 }
515 #endif
516 BFT_FREE(_var1_c);
517 BFT_FREE(_var2_c);
518
519 /* In case of single value, apply to all */
520
521 if (_n_vals == 1) {
522 cs_real_t val_const = val[0];
523 for (cs_lnum_t ii = 0; ii < n_vals; ii++)
524 val[ii] = val_const;
525 }
526 }
527
528 /*----------------------------------------------------------------------------*/
529 /*!
530 * \brief Compute properties with Freesteam in a defined thermal plane.
531 *
532 * \param[in] thermo_plane thermodynamic plane
533 * \param[in] property property queried
534 * \param[in] n_vals number of values
535 * \param[in] var1 values on first plane axis
536 * \param[in] var2 values on second plane axis
537 * \param[out] val resulting property values
538 */
539 /*----------------------------------------------------------------------------*/
540
541 void
cs_phys_prop_freesteam(cs_phys_prop_thermo_plane_type_t thermo_plane,cs_phys_prop_type_t property,const cs_lnum_t n_vals,const cs_real_t var1[],const cs_real_t var2[],cs_real_t val[])542 cs_phys_prop_freesteam(cs_phys_prop_thermo_plane_type_t thermo_plane,
543 cs_phys_prop_type_t property,
544 const cs_lnum_t n_vals,
545 const cs_real_t var1[],
546 const cs_real_t var2[],
547 cs_real_t val[])
548 {
549 #if defined(HAVE_FREESTEAM)
550 if (thermo_plane == CS_PHYS_PROP_PLANE_PH) {
551 for (cs_lnum_t i = 0; i < n_vals; i++) {
552 SteamState S0 = freesteam_set_ph(var1[i], var2[i]);
553 switch (property) {
554 case CS_PHYS_PROP_PRESSURE:
555 bft_error(__FILE__, __LINE__, 0,
556 _("bad choice: you choose to work in the %s plane."), "ph");
557 break;
558 case CS_PHYS_PROP_TEMPERATURE:
559 val[i] = freesteam_T(S0);
560 break;
561 case CS_PHYS_PROP_ENTHALPY:
562 bft_error(__FILE__, __LINE__, 0,
563 _("bad choice: you choose to work in the %s plane."), "ph");
564 break;
565 case CS_PHYS_PROP_ENTROPY:
566 val[i] = freesteam_s(S0);
567 break;
568 case CS_PHYS_PROP_ISOBARIC_HEAT_CAPACITY:
569 val[i] = freesteam_cp(S0);
570 break;
571 case CS_PHYS_PROP_ISOCHORIC_HEAT_CAPACITY:
572 val[i] = freesteam_cv(S0);
573 break;
574 case CS_PHYS_PROP_DENSITY:
575 val[i] = freesteam_rho(S0);
576 break;
577 case CS_PHYS_PROP_INTERNAL_ENERGY:
578 val[i] = freesteam_u(S0);
579 break;
580 case CS_PHYS_PROP_THERMAL_CONDUCTIVITY:
581 val[i] = freesteam_k(S0);
582 break;
583 case CS_PHYS_PROP_DYNAMIC_VISCOSITY:
584 val[i] = freesteam_mu(S0);
585 break;
586 case CS_PHYS_PROP_SPEED_OF_SOUND:
587 val[i] = freesteam_w(S0);
588 break;
589 }
590 }
591 }
592 else if (thermo_plane == CS_PHYS_PROP_PLANE_PT) {
593 for (cs_lnum_t i = 0; i < n_vals; i++) {
594 SteamState S0 = freesteam_set_pT(var1[i], var2[i]);
595 switch (property) {
596 case CS_PHYS_PROP_PRESSURE:
597 bft_error(__FILE__, __LINE__, 0,
598 _("bad choice: you choose to work in the %s plane."), "pT");
599 break;
600 case CS_PHYS_PROP_TEMPERATURE:
601 bft_error(__FILE__, __LINE__, 0,
602 _("bad choice: you choose to work in the %s plane."), "pT");
603 break;
604 case CS_PHYS_PROP_ENTHALPY:
605 val[i] = freesteam_h(S0);
606 break;
607 case CS_PHYS_PROP_ENTROPY:
608 val[i] = freesteam_s(S0);
609 break;
610 case CS_PHYS_PROP_ISOBARIC_HEAT_CAPACITY:
611 val[i] = freesteam_cp(S0);
612 break;
613 case CS_PHYS_PROP_ISOCHORIC_HEAT_CAPACITY:
614 val[i] = freesteam_cv(S0);
615 break;
616 case CS_PHYS_PROP_DENSITY:
617 val[i] = freesteam_rho(S0);
618 break;
619 case CS_PHYS_PROP_INTERNAL_ENERGY:
620 val[i] = freesteam_u(S0);
621 break;
622 case CS_PHYS_PROP_THERMAL_CONDUCTIVITY:
623 val[i] = freesteam_k(S0);
624 break;
625 case CS_PHYS_PROP_DYNAMIC_VISCOSITY:
626 val[i] = freesteam_mu(S0);
627 break;
628 case CS_PHYS_PROP_SPEED_OF_SOUND:
629 val[i] = freesteam_w(S0);
630 break;
631 }
632 }
633 }
634 else if (thermo_plane == CS_PHYS_PROP_PLANE_PS) {
635 for (cs_lnum_t i = 0; i < n_vals; i++) {
636 SteamState S0 = freesteam_set_ps(var1[i], var2[i]);
637 switch (property) {
638 case CS_PHYS_PROP_PRESSURE:
639 bft_error(__FILE__, __LINE__, 0,
640 _("bad choice: you choose to work in the %s plane."), "ps");
641 break;
642 case CS_PHYS_PROP_TEMPERATURE:
643 val[i] = freesteam_T(S0);
644 break;
645 case CS_PHYS_PROP_ENTHALPY:
646 val[i] = freesteam_h(S0);
647 break;
648 case CS_PHYS_PROP_ENTROPY:
649 bft_error(__FILE__, __LINE__, 0,
650 _("bad choice: you choose to work in the %s plane."), "ps");
651 break;
652 case CS_PHYS_PROP_ISOBARIC_HEAT_CAPACITY:
653 val[i] = freesteam_cp(S0);
654 break;
655 case CS_PHYS_PROP_ISOCHORIC_HEAT_CAPACITY:
656 val[i] = freesteam_cv(S0);
657 break;
658 case CS_PHYS_PROP_DENSITY:
659 val[i] = freesteam_rho(S0);
660 break;
661 case CS_PHYS_PROP_INTERNAL_ENERGY:
662 val[i] = freesteam_u(S0);
663 break;
664 case CS_PHYS_PROP_THERMAL_CONDUCTIVITY:
665 val[i] = freesteam_k(S0);
666 break;
667 case CS_PHYS_PROP_DYNAMIC_VISCOSITY:
668 val[i] = freesteam_mu(S0);
669 break;
670 case CS_PHYS_PROP_SPEED_OF_SOUND:
671 val[i] = freesteam_w(S0);
672 break;
673 }
674 }
675 }
676 else if (thermo_plane == CS_PHYS_PROP_PLANE_PU) {
677 for (cs_lnum_t i = 0; i < n_vals; i++) {
678 SteamState S0 = freesteam_set_pu(var1[i], var2[i]);
679 switch (property) {
680 case CS_PHYS_PROP_PRESSURE:
681 bft_error(__FILE__, __LINE__, 0,
682 _("bad choice: you choose to work in the %s plane."), "pu");
683 break;
684 case CS_PHYS_PROP_TEMPERATURE:
685 val[i] = freesteam_T(S0);
686 break;
687 case CS_PHYS_PROP_ENTHALPY:
688 val[i] = freesteam_h(S0);
689 break;
690 case CS_PHYS_PROP_ENTROPY:
691 val[i] = freesteam_s(S0);
692 break;
693 case CS_PHYS_PROP_ISOBARIC_HEAT_CAPACITY:
694 val[i] = freesteam_cp(S0);
695 break;
696 case CS_PHYS_PROP_ISOCHORIC_HEAT_CAPACITY:
697 val[i] = freesteam_cv(S0);
698 break;
699 case CS_PHYS_PROP_DENSITY:
700 val[i] = freesteam_rho(S0);
701 break;
702 case CS_PHYS_PROP_INTERNAL_ENERGY:
703 bft_error(__FILE__, __LINE__, 0,
704 _("bad choice: you choose to work in the %s plane."), "pu");
705 break;
706 case CS_PHYS_PROP_THERMAL_CONDUCTIVITY:
707 val[i] = freesteam_k(S0);
708 break;
709 case CS_PHYS_PROP_DYNAMIC_VISCOSITY:
710 val[i] = freesteam_mu(S0);
711 break;
712 case CS_PHYS_PROP_SPEED_OF_SOUND:
713 val[i] = freesteam_w(S0);
714 break;
715 }
716 }
717 }
718 else if (thermo_plane == CS_PHYS_PROP_PLANE_PV) {
719 for (cs_lnum_t i = 0; i < n_vals; i++) {
720 SteamState S0 = freesteam_set_pv(var1[i], var2[i]);
721 switch (property) {
722 case CS_PHYS_PROP_PRESSURE:
723 bft_error(__FILE__, __LINE__, 0,
724 _("bad choice: you choose to work in the %s plane."), "pv");
725 break;
726 case CS_PHYS_PROP_TEMPERATURE:
727 val[i] = freesteam_T(S0);
728 break;
729 case CS_PHYS_PROP_ENTHALPY:
730 val[i] = freesteam_h(S0);
731 break;
732 case CS_PHYS_PROP_ENTROPY:
733 val[i] = freesteam_s(S0);
734 break;
735 case CS_PHYS_PROP_ISOBARIC_HEAT_CAPACITY:
736 val[i] = freesteam_cp(S0);
737 break;
738 case CS_PHYS_PROP_ISOCHORIC_HEAT_CAPACITY:
739 val[i] = freesteam_cv(S0);
740 break;
741 case CS_PHYS_PROP_DENSITY:
742 val[i] = freesteam_rho(S0);
743 break;
744 case CS_PHYS_PROP_INTERNAL_ENERGY:
745 val[i] = freesteam_u(S0);
746 break;
747 case CS_PHYS_PROP_THERMAL_CONDUCTIVITY:
748 val[i] = freesteam_k(S0);
749 break;
750 case CS_PHYS_PROP_DYNAMIC_VISCOSITY:
751 val[i] = freesteam_mu(S0);
752 break;
753 case CS_PHYS_PROP_SPEED_OF_SOUND:
754 val[i] = freesteam_w(S0);
755 break;
756 }
757 }
758 }
759 else if (thermo_plane == CS_PHYS_PROP_PLANE_TS) {
760 for (cs_lnum_t i = 0; i < n_vals; i++) {
761 SteamState S0 = freesteam_set_Ts(var1[i], var2[i]);
762 switch (property) {
763 case CS_PHYS_PROP_PRESSURE:
764 val[i] = freesteam_p(S0);
765 break;
766 case CS_PHYS_PROP_TEMPERATURE:
767 bft_error(__FILE__, __LINE__, 0,
768 _("bad choice: you choose to work in the %s plane."), "Ts");
769 break;
770 case CS_PHYS_PROP_ENTHALPY:
771 val[i] = freesteam_h(S0);
772 break;
773 case CS_PHYS_PROP_ENTROPY:
774 bft_error(__FILE__, __LINE__, 0,
775 _("bad choice: you choose to work in the %s plane."), "Ts");
776 break;
777 case CS_PHYS_PROP_ISOBARIC_HEAT_CAPACITY:
778 val[i] = freesteam_cp(S0);
779 break;
780 case CS_PHYS_PROP_ISOCHORIC_HEAT_CAPACITY:
781 val[i] = freesteam_cv(S0);
782 break;
783 case CS_PHYS_PROP_DENSITY:
784 val[i] = freesteam_rho(S0);
785 break;
786 case CS_PHYS_PROP_INTERNAL_ENERGY:
787 val[i] = freesteam_u(S0);
788 break;
789 case CS_PHYS_PROP_THERMAL_CONDUCTIVITY:
790 val[i] = freesteam_k(S0);
791 break;
792 case CS_PHYS_PROP_DYNAMIC_VISCOSITY:
793 val[i] = freesteam_mu(S0);
794 break;
795 case CS_PHYS_PROP_SPEED_OF_SOUND:
796 val[i] = freesteam_w(S0);
797 break;
798 }
799 }
800 }
801 else if (thermo_plane == CS_PHYS_PROP_PLANE_TX) {
802 for (cs_lnum_t i = 0; i < n_vals; i++) {
803 SteamState S0 = freesteam_set_Tx(var1[i], var2[i]);
804 switch (property) {
805 case CS_PHYS_PROP_PRESSURE:
806 val[i] = freesteam_p(S0);
807 break;
808 case CS_PHYS_PROP_TEMPERATURE:
809 bft_error(__FILE__, __LINE__, 0,
810 _("bad choice: you choose to work in the %s plane."), "Tx");
811 break;
812 case CS_PHYS_PROP_ENTHALPY:
813 val[i] = freesteam_h(S0);
814 break;
815 case CS_PHYS_PROP_ENTROPY:
816 val[i] = freesteam_s(S0);
817 break;
818 case CS_PHYS_PROP_ISOBARIC_HEAT_CAPACITY:
819 val[i] = freesteam_cp(S0);
820 break;
821 case CS_PHYS_PROP_ISOCHORIC_HEAT_CAPACITY:
822 val[i] = freesteam_cv(S0);
823 break;
824 case CS_PHYS_PROP_DENSITY:
825 val[i] = freesteam_rho(S0);
826 break;
827 case CS_PHYS_PROP_INTERNAL_ENERGY:
828 val[i] = freesteam_u(S0);
829 break;
830 case CS_PHYS_PROP_THERMAL_CONDUCTIVITY:
831 val[i] = freesteam_k(S0);
832 break;
833 case CS_PHYS_PROP_DYNAMIC_VISCOSITY:
834 val[i] = freesteam_mu(S0);
835 break;
836 case CS_PHYS_PROP_SPEED_OF_SOUND:
837 val[i] = freesteam_w(S0);
838 break;
839 }
840 }
841 }
842 #else
843 CS_UNUSED(thermo_plane);
844 CS_UNUSED(property);
845 CS_UNUSED(n_vals);
846 CS_UNUSED(var1);
847 CS_UNUSED(var2);
848 CS_UNUSED(val);
849
850 bft_error(__FILE__, __LINE__, 0,
851 _("Freesteam support not available in this build."));
852 #endif
853 }
854
855
856 /*----------------------------------------------------------------------------*/
857 /*!
858 * \brief Get reference value of a physical property
859 *
860 * \param[in] name property name
861 *
862 * \return reference value (cs_real_t)
863 */
864 /*----------------------------------------------------------------------------*/
865
866 cs_real_t
cs_physical_property_get_ref_value(const char * name)867 cs_physical_property_get_ref_value(const char *name)
868 {
869
870 const cs_property_t *pty = cs_property_by_name(name);
871 if (pty == NULL)
872 bft_error(__FILE__, __LINE__, 0,
873 _("Error: property '%s' does not exist\n"),
874 name);
875
876 return pty->ref_value;
877
878 }
879
880
881 /*----------------------------------------------------------------------------*/
882 /*!
883 * \brief Set reference value for a physical property
884 *
885 * \param[in] name property name
886 * \param[in] val new value to set
887 */
888 /*----------------------------------------------------------------------------*/
889
890 void
cs_physical_property_set_ref_value(const char * name,const cs_real_t val)891 cs_physical_property_set_ref_value(const char *name,
892 const cs_real_t val)
893 {
894 cs_property_t *pty = cs_property_by_name(name);
895 if (pty == NULL)
896 bft_error(__FILE__, __LINE__, 0,
897 _("Error: property '%s' does not exist\n"),
898 name);
899
900 cs_property_set_reference_value(pty, val);
901 }
902
903 /*----------------------------------------------------------------------------*/
904 /*!
905 * \brief Get property reference values for a given zone
906 *
907 * \param[in] name property name
908 * \param[in] zname zone name
909 * \param[in] retval array of values to return
910 */
911 /*----------------------------------------------------------------------------*/
912
913 void
cs_physical_property_get_zone_values(const char * name,const char * zname,cs_real_t retval[])914 cs_physical_property_get_zone_values(const char *name,
915 const char *zname,
916 cs_real_t retval[])
917 {
918 const cs_property_t *pty = cs_property_by_name(name);
919 if (pty == NULL)
920 bft_error(__FILE__, __LINE__, 0,
921 _("Error: property '%s' does not exist\n"),
922 name);
923
924 cs_xdef_t *def = _get_property_def_on_zone(pty, zname);
925 if (def == NULL)
926 bft_error(__FILE__, __LINE__, 0,
927 _("Error: property '%s' does not have a definition for "
928 "zone '%s'\n"),
929 name, zname);
930
931 /* Sanity check */
932 assert(def->type == CS_XDEF_BY_VALUE);
933
934 if (pty->type & CS_PROPERTY_ISO) {
935 const cs_real_t *_context = (cs_real_t *)def->context;
936 retval[0] = _context[0];
937
938 } else if (pty->type & CS_PROPERTY_ORTHO) {
939 const cs_real_t *_context = (cs_real_t *)def->context;
940 for (int j = 0; j < 3; j++)
941 retval[j] = _context[j];
942
943 } else if (pty->type & CS_PROPERTY_ANISO_SYM) {
944 const cs_real_t *_context = (cs_real_t *)def->context;
945 for (int j = 0; j < 6; j++)
946 retval[j] = _context[j];
947
948 } else if (pty->type & CS_PROPERTY_ANISO) {
949 const cs_real_3_t *_context = (cs_real_3_t *)def->context;
950 for (int j = 0; j < 3; j++)
951 for (int k = 0; k < 3; k++)
952 retval[3*j + k] = _context[j][k];
953 }
954 }
955
956 /*----------------------------------------------------------------------------*/
957 /*!
958 * \brief Create a physical property
959 *
960 * \param[in] name property name
961 * \param[in] dim property dimension
962 * \param[in] refval reference value
963 */
964 /*----------------------------------------------------------------------------*/
965
966 void
cs_physical_property_create(const char * name,const int dim,const cs_real_t refval)967 cs_physical_property_create(const char *name,
968 const int dim,
969 const cs_real_t refval)
970 {
971 _physical_property_create(name, dim, refval);
972 }
973
974 /*----------------------------------------------------------------------------*/
975 /*!
976 * \brief Add a property definition on a given zone using a single value
977 *
978 * \param[in] name property name
979 * \param[in] zname zone name
980 * \param[in] dim property dimension
981 * \param[in] val reference value for the zone
982 */
983 /*----------------------------------------------------------------------------*/
984
985 void
cs_physical_property_define_from_value(const char * name,const char * zname,const int dim,const cs_real_t val)986 cs_physical_property_define_from_value(const char *name,
987 const char *zname,
988 const int dim,
989 const cs_real_t val)
990 {
991 cs_property_t *pty = cs_property_by_name(name);
992 if (pty == NULL)
993 pty = _physical_property_create(name, dim, 0.);
994
995 if (dim == 1) {
996 cs_property_def_iso_by_value(pty, zname, val);
997 }
998 else if (dim == 3) {
999 cs_real_t dvals[3] = {val, val, val};
1000 cs_property_def_ortho_by_value(pty, zname, dvals);
1001 } else if (dim == 6) {
1002 cs_real_t dvals[6] = {val, val, val, val, val, val};
1003 cs_property_def_aniso_sym_by_value(pty, zname, dvals);
1004 }
1005 else if (dim == 9) {
1006 cs_real_t dvals[3][3] = { {val, 0., 0.},
1007 {0., val, 0.},
1008 {0., 0., val} };
1009 cs_property_def_aniso_by_value(pty, zname, dvals);
1010 }
1011 }
1012
1013 /*----------------------------------------------------------------------------*/
1014 /*!
1015 * \brief Add a property multi-diemnsional definition on a given zone
1016 *
1017 * \param[in] name property name
1018 * \param[in] zname zone name
1019 * \param[in] dim property dimension (>1)
1020 * \param[in] vals array of values to set
1021 */
1022 /*----------------------------------------------------------------------------*/
1023
1024 void
cs_physical_property_define_from_values(const char * name,const char * zname,const int dim,cs_real_t vals[])1025 cs_physical_property_define_from_values(const char *name,
1026 const char *zname,
1027 const int dim,
1028 cs_real_t vals[])
1029 {
1030 assert(dim > 1 && vals != NULL);
1031
1032 cs_property_t *pty = cs_property_by_name(name);
1033
1034 if (pty == NULL)
1035 pty = _physical_property_create(name, dim, 0.);
1036
1037 if (dim == 3)
1038 cs_property_def_ortho_by_value(pty, zname, vals);
1039 else if (dim == 6)
1040 cs_property_def_aniso_sym_by_value(pty, zname, vals);
1041 else if (dim == 9) {
1042 cs_real_3_t *vals2use = (cs_real_3_t *)vals;
1043 cs_property_def_aniso_by_value(pty, zname, vals2use);
1044 }
1045 }
1046
1047 /*----------------------------------------------------------------------------*/
1048 /*!
1049 * \brief Add a property definition based on a cs_field_t.
1050 *
1051 * The field is created if needed
1052 *
1053 * \param[in] name property name
1054 * \param[in] type_flag field type flag
1055 * \param[in] location_id location id flag
1056 * \param[in] dim field dimension
1057 * \param[in] has_previous does the field has val_pre
1058 */
1059 /*----------------------------------------------------------------------------*/
1060
1061 void
cs_physical_property_define_from_field(const char * name,int type_flag,int location_id,int dim,bool has_previous)1062 cs_physical_property_define_from_field(const char *name,
1063 int type_flag,
1064 int location_id,
1065 int dim,
1066 bool has_previous)
1067 {
1068 cs_property_t *pty = cs_property_by_name(name);
1069 if (pty == NULL)
1070 pty = _physical_property_create(name, dim, 0.);
1071
1072 cs_field_t *f = cs_field_by_name_try(name);
1073 if (f == NULL)
1074 f = cs_field_create(name, type_flag, location_id, dim, has_previous);
1075
1076 cs_property_def_by_field(pty, f);
1077 }
1078
1079 /*----------------------------------------------------------------------------*/
1080 /*!
1081 * \brief Return id of field associated to property
1082 *
1083 * \param[in] name property name
1084 *
1085 * \return field id (int)
1086 */
1087 /*----------------------------------------------------------------------------*/
1088
1089 int
cs_physical_property_field_id_by_name(const char * name)1090 cs_physical_property_field_id_by_name(const char *name)
1091 {
1092 int retval = -1;
1093
1094 cs_field_t *f = cs_field_by_name_try(name);
1095
1096 if (f != NULL)
1097 retval = f->id;
1098
1099 return retval;
1100 }
1101
1102 /*----------------------------------------------------------------------------*/
1103 /*!
1104 * \brief Update reference values for a property on a given zone
1105 *
1106 * \param[in] name property name
1107 * \param[in] zname zone name
1108 * \param[in] vals array of values to set
1109 */
1110 /*----------------------------------------------------------------------------*/
1111
1112 void
cs_physical_property_update_zone_values(const char * name,const char * zname,const cs_real_t vals[])1113 cs_physical_property_update_zone_values(const char *name,
1114 const char *zname,
1115 const cs_real_t vals[])
1116 {
1117 cs_property_t *pty = cs_property_by_name(name);
1118
1119 cs_xdef_t *def = _get_property_def_on_zone(pty, zname);
1120
1121 _update_def_values(def, vals);
1122 }
1123
1124
1125 /*----------------------------------------------------------------------------*/
1126
1127 END_C_DECLS
1128