sgas.c (19224B)
1 /* Copyright (C) 2025, 2026 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2025, 2026 Université de Lorraine 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #define _POSIX_C_SOURCE 200112L /* strtok_r and wordexp */ 18 19 #include "sgas.h" 20 21 #include <astoria/atrtp.h> 22 #include <star/smsh.h> 23 #include <star/suvm.h> 24 25 #include <rsys/algorithm.h> 26 #include <rsys/clock_time.h> 27 #include <rsys/cstr.h> 28 #include <rsys/dynamic_array.h> 29 #include <rsys/logger.h> 30 #include <rsys/mem_allocator.h> 31 #include <rsys/ref_count.h> 32 #include <rsys/text_reader.h> 33 34 #include <string.h> 35 #include <wordexp.h> 36 37 /* Helper macros for logging */ 38 #define LOG__(Gas, Lvl, Type, ...) { \ 39 if ((Gas)->verbose >= (Lvl)) \ 40 logger_print((Gas)->logger, Type, __VA_ARGS__); \ 41 } (void)0 42 #define ERROR(Gas, ...) LOG__(Gas, 1, LOG_ERROR, "error: "__VA_ARGS__) 43 #define WARN(Gas, ...) LOG__(Gas, 2, LOG_WARNING, "warning: "__VA_ARGS__) 44 #define INFO(Gas, ...) LOG__(Gas, 3, LOG_OUTPUT, __VA_ARGS__) 45 46 struct therm_props { 47 struct atrtp* props; 48 struct atrtp_desc desc; 49 double time; /* [s] ??? Confirm with Yaniss */ 50 }; 51 #define THERM_PROPS_NULL__ {0} 52 static const struct therm_props THERM_PROPS_NULL = THERM_PROPS_NULL__; 53 54 static INLINE void 55 therm_props_init(struct mem_allocator* allocator, struct therm_props* props) 56 { 57 ASSERT(props); 58 (void)allocator; 59 *props = THERM_PROPS_NULL; 60 } 61 62 static INLINE void 63 therm_props_release(struct therm_props* props) 64 { 65 ASSERT(props); 66 if(props->props) ATRTP(ref_put(props->props)); 67 } 68 69 static INLINE res_T 70 therm_props_copy(struct therm_props* dst, const struct therm_props* src) 71 { 72 ASSERT(src && dst); 73 if(src == dst) return RES_OK; 74 75 if(src->props) ATRTP(ref_get(src->props)); 76 if(dst->props) ATRTP(ref_put(dst->props)); 77 dst->props = src->props; 78 dst->desc = src->desc; 79 dst->time = src->time; 80 return RES_OK; 81 } 82 83 static INLINE res_T 84 therm_props_copy_and_release(struct therm_props* dst, struct therm_props* src) 85 { 86 ASSERT(src && dst); 87 88 if(dst->props) ATRTP(ref_put(dst->props)); 89 90 dst->props = src->props; 91 dst->desc = src->desc; 92 dst->time = src->time; 93 94 *src = THERM_PROPS_NULL; 95 96 return RES_OK; 97 } 98 99 /* Define the dynamic array of thermodynamic properties */ 100 #define DARRAY_NAME therm_props 101 #define DARRAY_DATA struct therm_props 102 #define DARRAY_FUNCTOR_INIT therm_props_init 103 #define DARRAY_FUNCTOR_RELEASE therm_props_release 104 #define DARRAY_FUNCTOR_COPY therm_props_copy 105 #define DARRAY_FUNCTOR_COPY_AND_RELEASE therm_props_copy_and_release 106 #include <rsys/dynamic_array.h> 107 108 struct sgas { 109 /* Structured geometry */ 110 struct suvm_volume* mesh; 111 struct suvm_mesh_desc mesh_desc; 112 113 /* Time varying thermodynamic properties */ 114 struct darray_therm_props therm_props; 115 116 struct mem_allocator* allocator; 117 struct logger* logger; 118 int verbose; /* Verbosity level */ 119 ref_T ref; 120 }; 121 122 /******************************************************************************* 123 * Helper functions 124 ******************************************************************************/ 125 static INLINE res_T 126 check_sgas_file(const struct sgas_file* file) 127 { 128 return (!file || (!file->name && !file->stream)) ? RES_BAD_ARG : RES_OK; 129 } 130 131 static INLINE res_T 132 check_sgas_create_args(const struct sgas_create_args* args) 133 { 134 res_T res = RES_OK; 135 136 if(!args) return RES_BAD_ARG; 137 if((res = check_sgas_file(&args->mesh)) != RES_OK) return res; 138 if((res = check_sgas_file(&args->therm_props_list)) != RES_OK) return res; 139 140 return RES_OK; 141 } 142 143 static INLINE res_T 144 check_therm_props 145 (struct sgas* gas, 146 /* Thermo Properties to be checked against the gas mesh 147 * and previously recorded ones properties */ 148 const struct therm_props* props, 149 /* Name of the file listing the thermodynamic properties and line in which 150 * the thermodynamic properties are defined */ 151 const char* list_filename, 152 const size_t list_line, 153 /* Name of the file from which the thermodynamic properties were loaded */ 154 const char* atrtp_filename) 155 { 156 size_t n = 0; 157 res_T res = RES_OK; 158 ASSERT(gas && props && list_filename && atrtp_filename); 159 160 /* Check that thermodynamic properties are defined per tetrahedron. 161 * 162 * Note that the atrtp descriptor assumes that thermodynamic data is defined 163 * per node, not functionally, but using the name of the “node” as the a 164 * priori data carrier. But in practice, the only constraint is that the 165 * data indexing does not exceed the total number of data points. Thus, 166 * “nodes” can in practice be “cells.” 167 * 168 * TODO In any case, this is confusing, so the ATRTP API needs to be updated 169 * so that the data indexing appears to be independent of the data with 170 * which it is associated. */ 171 if(gas->mesh_desc.nprimitives != props->desc.nnodes) { 172 ERROR(gas, 173 "%s: the number of thermodynamic property sets does not correspond " 174 "to the number of tetrahedra in the gas mesh.\n", 175 atrtp_filename); 176 res = RES_BAD_ARG; 177 goto error; 178 } 179 180 n = darray_therm_props_size_get(&gas->therm_props); 181 if(n) { 182 const struct therm_props* props_prev = NULL; 183 184 /* Check that the the therm properties are listed in chronological order */ 185 props_prev = darray_therm_props_cdata_get(&gas->therm_props)+n-1; 186 if(props->time < props_prev->time) { 187 ERROR(gas, 188 "%s:%lu: the thermodynamic properties must be listed " 189 "in chronological order\n", list_filename, (unsigned long)list_line); 190 res = RES_BAD_ARG; 191 goto error; 192 } 193 } 194 195 exit: 196 return res; 197 error: 198 goto exit; 199 } 200 201 static INLINE int 202 cmp_therm_props_time(const void* a, const void* b) 203 { 204 const struct therm_props* props = NULL; 205 double time = 0; 206 ASSERT(a && b); 207 208 time = *(double*)(a); 209 props = b; 210 211 if(time < props->time) { 212 return -1; 213 } else if(time > props->time) { 214 return +1; 215 } else { 216 return 0; 217 } 218 } 219 220 static const struct therm_props* 221 fetch_therm_props(const struct sgas* gas, const double time) 222 { 223 const struct therm_props* find = NULL; 224 const struct therm_props* props = NULL; 225 const struct therm_props* prop_list = NULL; 226 size_t nprops = 0; 227 228 ASSERT(gas); 229 230 prop_list = darray_therm_props_cdata_get(&gas->therm_props); 231 nprops = darray_therm_props_size_get(&gas->therm_props); 232 233 /* Search for the first set of thermodynamic properties whose time is greater 234 * than or equal to the input time */ 235 find = search_lower_bound 236 (&time, prop_list, nprops, sizeof(*prop_list), cmp_therm_props_time); 237 238 if(find == prop_list) { 239 /* The input time is earlier than that of the first set of thermodynamic 240 * properties. Query this first entry */ 241 props = find; 242 243 } else if(!find) { 244 /* The input time is later than that of the last set of thermodynamic 245 * properties. Query this last entry. */ 246 props = prop_list + nprops-1; 247 248 } else { 249 /* The input time lies between two sets of thermodynamic properties. Query 250 * the one that is closest. */ 251 ASSERT(find->time >= time && find != prop_list); 252 props = (find[0].time - time) < (time - find[-1].time) ? find : find-1; 253 } 254 255 return props; 256 } 257 258 static void 259 tetra_mesh_get_indices(const size_t itetra, size_t out_ids[4], void* ctx) 260 { 261 const struct smsh_desc* desc = ctx; 262 const uint64_t* ids = NULL; 263 264 /* Pre-conditions */ 265 ASSERT(desc && out_ids && itetra < desc->ncells && desc->dcell == 4); 266 267 ids = smsh_desc_get_cell(desc, itetra); 268 out_ids[0] = ids[0]; 269 out_ids[1] = ids[1]; 270 out_ids[2] = ids[2]; 271 out_ids[3] = ids[3]; 272 } 273 274 static void 275 tetra_mesh_get_position(const size_t ivert, double out_pos[3], void* ctx) 276 { 277 const struct smsh_desc* desc = ctx; 278 const double* pos = NULL; 279 280 /* Pre-conditions */ 281 ASSERT(desc && out_pos && ivert < desc->nnodes && desc->dnode == 3); 282 283 pos = smsh_desc_get_node(desc, ivert); 284 out_pos[0] = pos[0]; 285 out_pos[1] = pos[1]; 286 out_pos[2] = pos[2]; 287 } 288 289 static res_T 290 build_accel_struct 291 (struct sgas* gas, 292 const struct smsh* mesh) /* Unstructured mesh */ 293 { 294 struct suvm_tetrahedral_mesh_args mesh_args = SUVM_TETRAHEDRAL_MESH_ARGS_NULL; 295 struct suvm_device* suvm = NULL; 296 struct smsh_desc mesh_desc = SMSH_DESC_NULL; 297 res_T res = RES_OK; 298 ASSERT(gas && mesh); /* medium */ 299 300 res = suvm_device_create 301 (gas->logger, 302 gas->allocator, 303 gas->verbose, 304 &suvm); 305 if(res != RES_OK) goto error; 306 307 res = smsh_get_desc(mesh, &mesh_desc); 308 if(res != RES_OK) goto error; 309 310 mesh_args.ntetrahedra = mesh_desc.ncells; 311 mesh_args.nvertices = mesh_desc.nnodes; 312 mesh_args.get_indices = tetra_mesh_get_indices; 313 mesh_args.get_position = tetra_mesh_get_position; 314 mesh_args.context = &mesh_desc; 315 res = suvm_tetrahedral_mesh_create(suvm, &mesh_args, &gas->mesh); 316 if(res != RES_OK) goto error; 317 318 res = suvm_volume_get_mesh_desc(gas->mesh, &gas->mesh_desc); 319 if(res != RES_OK) goto error; 320 321 exit: 322 if(suvm) SUVM(device_ref_put(suvm)); 323 return res; 324 error: 325 ERROR(gas, "Gas mesh structuring error -- %s\n", 326 res_to_cstr(res)); 327 goto exit; 328 } 329 330 static res_T 331 load_mesh 332 (struct sgas* gas, 333 const struct sgas_create_args* args, 334 struct smsh** out_mesh) 335 { 336 #define MESH_NAME (args->mesh.name ? args->mesh.name : "stream") 337 338 struct smsh_create_args create_args = SMSH_CREATE_ARGS_DEFAULT; 339 struct smsh_desc mesh_desc = SMSH_DESC_NULL; 340 struct smsh* mesh = NULL; 341 res_T res = RES_OK; 342 343 ASSERT(gas && args && out_mesh); /* Pre-conditions */ 344 345 create_args.allocator = gas->allocator; 346 create_args.logger = gas->logger; 347 create_args.verbose = gas->verbose; 348 if((res = smsh_create(&create_args, &mesh)) != RES_OK) goto error; 349 350 351 if(!args->mesh.stream) { 352 /* Load from file */ 353 struct smsh_load_args load_args = SMSH_LOAD_ARGS_NULL; 354 load_args.path = MESH_NAME; 355 load_args.memory_mapping = 1; 356 if((res = smsh_load(mesh, &load_args)) != RES_OK) goto error; 357 358 } else { 359 /* Load from stream */ 360 struct smsh_load_stream_args load_args = SMSH_LOAD_STREAM_ARGS_NULL; 361 load_args.stream = args->mesh.stream; 362 load_args.name = MESH_NAME; 363 load_args.memory_mapping = 1; 364 if((res = smsh_load_stream(mesh, &load_args)) != RES_OK) goto error; 365 } 366 367 res = smsh_get_desc(mesh, &mesh_desc); 368 if(res != RES_OK) goto error; 369 370 if(mesh_desc.dcell != 4) { 371 ERROR(gas, "%s: it is not a tetrahedral mesh as expected\n", MESH_NAME); 372 res = RES_BAD_ARG; 373 goto error; 374 } 375 376 exit: 377 *out_mesh = mesh; 378 return res; 379 error: 380 ERROR(gas, "Error loading gas mesh \"%s\" -- %s\n", 381 MESH_NAME, res_to_cstr(res)); 382 if(mesh) { SMSH(ref_put(mesh)); mesh = NULL; } 383 goto exit; 384 385 #undef MESH_NAME 386 } 387 388 static res_T 389 setup_mesh(struct sgas* gas, const struct sgas_create_args* args) 390 { 391 struct smsh* mesh = NULL; 392 res_T res = RES_OK; 393 394 ASSERT(gas && args); /* Pre-conditions */ 395 396 res = load_mesh(gas, args, &mesh); 397 if(res != RES_OK) goto error; 398 399 res = build_accel_struct(gas, mesh); 400 if(res != RES_OK) goto error; 401 402 exit: 403 if(mesh) SMSH(ref_put(mesh)); 404 return res; 405 error: 406 goto exit; 407 } 408 409 static res_T 410 parse_therm_props(struct sgas* gas, struct txtrdr* txtrdr) 411 { 412 /* Expansion */ 413 wordexp_t wexp; 414 int wexp_is_allocated = 0; 415 int err = 0; 416 417 /* Input file */ 418 const char* fname = NULL; /* Name of the read file */ 419 size_t lnum = 0; /* Current parsed line */ 420 421 /* Parsed line */ 422 char* line = NULL; 423 char* tk = NULL; 424 char* tk_ctx = NULL; 425 426 /* Miscellaneous */ 427 struct therm_props therm_props; 428 res_T res = RES_OK; 429 ASSERT(gas && txtrdr); 430 431 therm_props_init(gas->allocator, &therm_props); 432 433 /* Recover the file name that lists thermodynamic properties that vary over 434 * time, and the line (content and number) currently analyzed in this file, 435 * which defines the thermodynamic properties file currently analyzed. */ 436 fname = txtrdr_get_name(txtrdr); 437 lnum = txtrdr_get_line_num(txtrdr); 438 line = txtrdr_get_line(txtrdr); 439 ASSERT(line); 440 441 /* Parse the time */ 442 tk = strtok_r(line, " \t", &tk_ctx); 443 ASSERT(tk); /* The line is not empty and should contain at least one string */ 444 res = cstr_to_double(tk, &therm_props.time); 445 if(res != RES_OK) { 446 ERROR(gas, "%s:%lu: invalid time \"%s\"\n", fname, lnum, tk); 447 goto error; 448 } 449 450 /* Parse the file name of the thermodynamic properties */ 451 tk = strtok_r(NULL, "", &tk_ctx); /* Get the rest of the line */ 452 if(!tk) { 453 ERROR(gas, "%s:%lu: the thermodynamic properties are missing\n", 454 fname, lnum); 455 res = RES_BAD_ARG; 456 goto error; 457 } 458 459 /* Expands variables in file name, i.e. substitute their value */ 460 err = wordexp(tk, &wexp, 0/*flags*/); 461 if(err) { 462 ERROR(gas, "%s:%lu: unable to expand the filename string\n", 463 fname, lnum); 464 res = RES_BAD_ARG; 465 goto error; 466 } 467 wexp_is_allocated = 1; /* Save a word expansion has been allocated */ 468 ASSERT(wexp.we_wordc != 0); 469 470 /* Check that the name of the path to thermodynamic properties has no space, 471 * or, more precisely, no unprotected space, i.e. not preceded by a backslash. 472 * Otherwise, warn the user that the file name is certainly poorly formed. */ 473 if(wexp.we_wordc > 1) { 474 WARN(gas, 475 "%s:%lu: unexpected space in the path to the thermodynamic properties\n", 476 fname, lnum); 477 } 478 479 /* Here we go! 480 * 1 - Create the atrp device */ 481 res = atrtp_create 482 (gas->logger, 483 gas->allocator, 484 gas->verbose, 485 &therm_props.props); 486 if(res != RES_OK) goto error; 487 488 /* 2 - Load the thermodynamic properties */ 489 res = atrtp_load(therm_props.props, wexp.we_wordv[0]); 490 if(res != RES_OK) goto error; 491 492 /* 3 - Retrieve their descriptor */ 493 res = atrtp_get_desc(therm_props.props, &therm_props.desc); 494 if(res != RES_OK) goto error; 495 496 /* 4 - Check them against the already loaded data */ 497 res = check_therm_props(gas, &therm_props, fname, lnum, tk/*atrtp filename*/); 498 if(res != RES_OK) goto error; 499 500 /* 5 - Save them in the list of available thermo dynamic properties */ 501 res = darray_therm_props_push_back(&gas->therm_props, &therm_props); 502 if(res != RES_OK) { 503 ERROR(gas, 504 "%s:%lu: error when saving thermodynamic properties -- %s\n", 505 fname, lnum, res_to_cstr(res)); 506 goto error; 507 } 508 509 exit: 510 if(wexp_is_allocated) wordfree(&wexp); 511 therm_props_release(&therm_props); 512 return res; 513 error: 514 goto exit; 515 } 516 517 static res_T 518 setup_therm_props 519 (struct sgas* gas, 520 const struct sgas_create_args* args) 521 { 522 struct txtrdr* txtrdr; 523 res_T res = RES_OK; 524 525 ASSERT(gas && args); /* Pre-conditions */ 526 527 #define PROPS_NAME \ 528 (args->therm_props_list.name ? args->therm_props_list.name : "stream") 529 530 if(!args->therm_props_list.stream) { 531 /* Read the list of thermodynamic properties from an input file */ 532 res = txtrdr_file(gas->allocator, PROPS_NAME, '#', &txtrdr); 533 534 } else { 535 /* Read the list of thermodynamic properties from a stream */ 536 res = txtrdr_stream 537 (gas->allocator, args->therm_props_list.stream, PROPS_NAME, '#', &txtrdr); 538 } 539 if(res != RES_OK) { 540 ERROR(gas, "Cannot create \"%s\" reader -- %s\n", 541 PROPS_NAME, res_to_cstr(res)); 542 goto error; 543 } 544 545 for(;;) { 546 res = txtrdr_read_line(txtrdr); 547 if(res != RES_OK) { 548 ERROR(gas, "%s:%lu error reading the line -- %s\n", 549 txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr), 550 res_to_cstr(res)); 551 goto error; 552 } 553 554 if(!txtrdr_get_cline(txtrdr)) break; /* No more line to parse */ 555 556 res = parse_therm_props(gas, txtrdr); 557 if(res != RES_OK) goto error; 558 } 559 560 if(darray_therm_props_size_get(&gas->therm_props) == 0) { 561 ERROR(gas, "%s: the list of thermodynamic properties is empty\n", 562 txtrdr_get_name(txtrdr)); 563 res = RES_BAD_ARG; 564 goto error; 565 } 566 567 #undef PROPS_NAME 568 569 exit: 570 if(txtrdr) txtrdr_ref_put(txtrdr); 571 return res; 572 error: 573 darray_therm_props_clear(&gas->therm_props); 574 goto exit; 575 } 576 577 static void 578 release_gas(ref_T* ref) 579 { 580 struct sgas* gas = CONTAINER_OF(ref, struct sgas, ref); 581 ASSERT(ref); 582 583 if(gas->mesh) SUVM(volume_ref_put(gas->mesh)); 584 darray_therm_props_release(&gas->therm_props); 585 MEM_RM(gas->allocator, gas); 586 } 587 588 /******************************************************************************* 589 * Local functions 590 ******************************************************************************/ 591 res_T 592 sgas_create 593 (const struct sgas_create_args* args, 594 struct sgas** out_gas) 595 { 596 char buf[512]; 597 struct time t0, t1, elapsed_time; 598 struct mem_allocator* allocator = NULL; 599 struct sgas* gas = NULL; 600 res_T res = RES_OK; 601 602 time_current(&t0); 603 604 if(!args || !out_gas) { res = RES_BAD_ARG; goto error; } 605 res = check_sgas_create_args(args); 606 if(res != RES_OK) goto error; 607 608 allocator = args->allocator ? args->allocator : &mem_default_allocator; 609 610 gas = MEM_CALLOC(allocator, 1, sizeof(*gas)); 611 if(!gas) { res = RES_MEM_ERR; goto error; } 612 ref_init(&gas->ref); 613 gas->allocator = allocator; 614 gas->logger = gas->logger ? gas->logger : LOGGER_DEFAULT; 615 gas->verbose = args->verbose; 616 darray_therm_props_init(gas->allocator, &gas->therm_props); 617 618 if((res = setup_mesh(gas, args)) != RES_OK) goto error; 619 if((res = setup_therm_props(gas, args)) != RES_OK) goto error; 620 time_current(&t1); 621 time_sub(&elapsed_time, &t1, &t0); 622 623 time_dump(&elapsed_time, TIME_ALL, NULL, buf, sizeof(buf)); 624 INFO(gas, "Load and structure gas data in %s\n", buf); 625 626 exit: 627 if(out_gas) *out_gas = gas; 628 return res; 629 error: 630 if(gas) { sgas_ref_put(gas); gas = NULL; } 631 goto exit; 632 } 633 634 res_T 635 sgas_ref_get(struct sgas* gas) 636 { 637 if(!gas) return RES_BAD_ARG; 638 ref_get(&gas->ref); 639 return RES_OK; 640 } 641 642 res_T 643 sgas_ref_put(struct sgas* gas) 644 { 645 if(!gas) return RES_BAD_ARG; 646 ref_put(&gas->ref, release_gas); 647 return RES_OK; 648 } 649 650 res_T 651 sgas_get_props 652 (const struct sgas* gas, 653 const double pos[3], 654 const double time, 655 struct sgas_props* gas_props) 656 { 657 struct suvm_primitive prim = SUVM_PRIMITIVE_NULL; 658 double bcoords[4] = {0,0,0,0}; 659 res_T res = RES_OK; 660 661 if(!gas || !pos || !gas_props) { res = RES_BAD_ARG; goto error; } 662 663 res = suvm_volume_at(gas->mesh, pos, &prim, bcoords); 664 if(res != RES_OK) goto error; 665 666 if(SUVM_PRIMITIVE_NONE(&prim)) { /* Is the position outside the gas? */ 667 *gas_props = SGAS_PROPS_NULL; 668 669 } else { 670 const double* raw_props = NULL; 671 const struct therm_props* therm_props = NULL; 672 673 therm_props = fetch_therm_props(gas, time); 674 675 /* Note that the function is missnamed. The properties are retrieved from 676 * the _cell_ not from the node. */ 677 raw_props = atrtp_desc_get_node_properties(&therm_props->desc, prim.iprim); 678 679 gas_props->pressure = raw_props[ATRTP_PRESSURE]; 680 gas_props->temperature = raw_props[ATRTP_TEMPERATURE]; 681 gas_props->xH2O = raw_props[ATRTP_XH20]; 682 gas_props->xCO2 = raw_props[ATRTP_XCO2]; 683 gas_props->xCO = raw_props[ATRTP_XCO]; 684 } 685 686 exit: 687 return res; 688 error: 689 goto exit; 690 }