sln_tree.c (19788B)
1 /* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2026 Université de Lorraine 3 * Copyright (C) 2022 Centre National de la Recherche Scientifique 4 * Copyright (C) 2022 Université Paul Sabatier 5 * 6 * This program is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation, either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 18 19 #include "sln.h" 20 #include "sln_device_c.h" 21 #include "sln_line.h" 22 #include "sln_tree_c.h" 23 24 #include <star/shtr.h> 25 26 #include <rsys/algorithm.h> 27 #include <rsys/cstr.h> 28 #include <rsys/math.h> 29 30 struct stream { 31 const char* name; 32 FILE* fp; 33 int intern_fp; /* Define if the stream was internally opened */ 34 }; 35 static const struct stream STREAM_NULL = {NULL, NULL, 0}; 36 37 /******************************************************************************* 38 * Helper functions 39 ******************************************************************************/ 40 /* Check that the sum of the molecular concentrations is equal to 1 */ 41 static res_T 42 check_molecule_concentration 43 (struct sln_device* sln, 44 const char* caller, 45 const struct sln_tree_create_args* args) 46 { 47 int i = 0; 48 double sum = 0; 49 ASSERT(sln && caller && args); 50 51 FOR_EACH(i, 0, SHTR_MAX_MOLECULE_COUNT) { 52 if(i == SHTR_MOLECULE_ID_NULL) continue; 53 sum += args->molecules[i].concentration; 54 } 55 56 /* The sum of molecular concentrations must be less than or equal to 1. It may 57 * be less than 1 if the remaining part of the mixture is (implicitly) defined 58 * as a radiatively inactive gas */ 59 if(sum > 1 && sum-1 > 1e-6) { 60 ERROR(sln, 61 "%s: the sum of molecule concentrations is greater than 1: %g\n", 62 caller, sum); 63 return RES_BAD_ARG; 64 } 65 66 return RES_OK; 67 } 68 69 /* Verify that the isotope abundance are valids */ 70 static res_T 71 check_molecule_isotope_abundances 72 (struct sln_device* sln, 73 const char* caller, 74 const struct sln_molecule* molecule) 75 { 76 int i = 0; 77 double sum = 0; 78 ASSERT(sln && caller && molecule); 79 80 /* The isotopic abundances are the default ones. Nothing to do */ 81 if(!molecule->non_default_isotope_abundances) return RES_OK; 82 83 /* The isotopic abundances are not the default ones. 84 * Verify that they are valid ... */ 85 FOR_EACH(i, 0, SLN_MAX_ISOTOPES_COUNT) { 86 if(molecule->isotopes[i].abundance < 0) { 87 const int isotope_id = i + 1; /* isotope id in [1, 10] */ 88 ERROR(sln, "%s: invalid abundance of isotopie %d of %s: %g.\n", 89 caller, isotope_id, shtr_molecule_cstr(i), 90 molecule->isotopes[i].abundance); 91 return RES_BAD_ARG; 92 } 93 94 sum += molecule->isotopes[i].abundance; 95 } 96 97 /* ... and that their sum equals 1 */ 98 if(!eq_eps(sum, 1, 1e-6)) { 99 ERROR(sln, "%s: the %s isotope abundances does not sum to 1: %g.\n", 100 caller, shtr_molecule_cstr(i), sum); 101 return RES_BAD_ARG; 102 } 103 104 return RES_OK; 105 } 106 107 static res_T 108 check_molecules 109 (struct sln_device* sln, 110 const char* caller, 111 const struct sln_tree_create_args* args) 112 { 113 char molecule_ok[SHTR_MAX_MOLECULE_COUNT] = {0}; 114 115 double concentrations_sum = 0; 116 size_t iline = 0; 117 size_t nlines = 0; 118 res_T res = RES_OK; 119 ASSERT(args->lines); 120 121 res = check_molecule_concentration(sln, caller, args); 122 if(res != RES_OK) return res; 123 124 /* Iterate over the lines to define which molecules has to be checked, i.e., 125 * the ones used in the mixture */ 126 SHTR(line_list_get_size(args->lines, &nlines)); 127 FOR_EACH(iline, 0, nlines) { 128 struct shtr_line line = SHTR_LINE_NULL; 129 const struct sln_molecule* molecule = NULL; 130 131 SHTR(line_list_at(args->lines, iline, &line)); 132 133 /* This molecule was already checked */ 134 if(molecule_ok[line.molecule_id]) continue; 135 136 molecule = args->molecules + line.molecule_id; 137 138 if(molecule->concentration == 0) { 139 /* A molecular concentration of zero is allowed, but may be a user error, 140 * as 0 is the default concentration in the tree creation arguments. 141 * Therefore, warn the user about this value so that they can determine 142 * whether or not it is an error on their part. */ 143 WARN(sln, "%s: the concentration of %s is zero\n.\n", 144 caller, shtr_molecule_cstr(line.molecule_id)); 145 146 } else if(molecule->concentration < 0) { 147 /* Concentration cannot be negative... */ 148 ERROR(sln, "%s: invalid %s concentration: %g.\n", 149 FUNC_NAME, shtr_molecule_cstr(line.molecule_id), 150 molecule->concentration); 151 return RES_BAD_ARG; 152 } 153 154 concentrations_sum += molecule->concentration; 155 156 if(molecule->cutoff <= 0) { 157 /* ... cutoff either */ 158 ERROR(sln, "%s: invalid %s cutoff: %g.\n", 159 caller, shtr_molecule_cstr(line.molecule_id), molecule->cutoff); 160 return RES_BAD_ARG; 161 } 162 163 res = check_molecule_isotope_abundances(sln, caller, molecule); 164 if(res != RES_OK) return res; 165 166 molecule_ok[line.molecule_id] = 1; 167 } 168 169 /* The sum of molecular concentrations must be less than or equal to 1. It may 170 * be less than 1 if the remaining part of the mixture is (implicitly) defined 171 * as a radiatively inactive gas */ 172 if(concentrations_sum > 1 && (concentrations_sum - 1) > 1e-6) { 173 ERROR(sln, 174 "%s: the sum of molecule concentrations is greater than 1: %g\n", 175 caller, concentrations_sum); 176 return RES_BAD_ARG; 177 } 178 179 return RES_OK; 180 } 181 182 static res_T 183 check_sln_tree_create_args 184 (struct sln_device* sln, 185 const char* caller, 186 const struct sln_tree_create_args* args) 187 { 188 res_T res = RES_OK; 189 ASSERT(sln && caller); 190 191 if(!args) return RES_BAD_ARG; 192 193 if(!args->metadata) { 194 ERROR(sln, "%s: the metadata is missing.\n", caller); 195 return RES_BAD_ARG; 196 } 197 198 if(!args->lines) { 199 ERROR(sln, "%s: the list of lines is missing.\n", caller); 200 return RES_BAD_ARG; 201 } 202 203 if(args->nvertices_hint == 0) { 204 ERROR(sln, 205 "%s: invalid hint on the number of vertices around the line center %lu.\n", 206 caller, (unsigned long)args->nvertices_hint); 207 return RES_BAD_ARG; 208 } 209 210 if(args->mesh_decimation_err < 0) { 211 ERROR(sln, "%s: invalid decimation error %g.\n", 212 caller, args->mesh_decimation_err); 213 return RES_BAD_ARG; 214 } 215 216 if((unsigned)args->mesh_type >= SLN_MESH_TYPES_COUNT__) { 217 ERROR(sln, "%s: invalid mesh type %d.\n", caller, args->mesh_type); 218 return RES_BAD_ARG; 219 } 220 221 if((unsigned)args->line_profile >= SLN_LINE_PROFILES_COUNT__) { 222 ERROR(sln, "%s: invalid line profile %d.\n", caller, args->line_profile); 223 return RES_BAD_ARG; 224 } 225 226 res = check_molecules(sln, caller, args); 227 if(res != RES_OK) return res; 228 229 return RES_OK; 230 } 231 232 static res_T 233 check_sln_tree_read_args 234 (struct sln_device* sln, 235 const char* caller, 236 const struct sln_tree_read_args* args) 237 { 238 if(!args) return RES_BAD_ARG; 239 240 if(!args->shtr) { 241 ERROR(sln, 242 "%s: the handle to the Star-HITRAN library is missing.\n", caller); 243 return RES_BAD_ARG; 244 } 245 246 if(!args->file && !args->filename) { 247 ERROR(sln, 248 "%s: the source file is missing. No file name or stream is provided.\n", 249 caller); 250 return RES_BAD_ARG; 251 } 252 253 return RES_OK; 254 } 255 256 static res_T 257 check_sln_tree_write_args 258 (struct sln_device* sln, 259 const char* caller, 260 const struct sln_tree_write_args* args) 261 { 262 if(!args) return RES_BAD_ARG; 263 264 if(!args->file && !args->filename) { 265 ERROR(sln, 266 "%s: the destination file is missing. " 267 "No file name or stream is provided.\n", 268 caller); 269 return RES_BAD_ARG; 270 } 271 272 return RES_OK; 273 } 274 static INLINE void 275 stream_release(struct stream* stream) 276 { 277 ASSERT(stream); 278 if(stream->intern_fp && stream->fp) CHK(fclose(stream->fp) == 0); 279 } 280 281 static res_T 282 stream_init 283 (struct sln_device* sln, 284 const char* caller, 285 const char* name, /* NULL <=> default stream name */ 286 FILE* fp, /* NULL <=> open file "name" */ 287 const char* mode, /* mode in fopen */ 288 struct stream* stream) 289 { 290 res_T res = RES_OK; 291 292 ASSERT(sln && caller && stream); 293 ASSERT(fp || (name && mode)); 294 295 *stream = STREAM_NULL; 296 297 if(fp) { 298 stream->intern_fp = 0; 299 stream->name = name ? name : "stream"; 300 stream->fp = fp; 301 302 } else { 303 stream->intern_fp = 1; 304 stream->name = name; 305 if(!(stream->fp = fopen(name, mode))) { 306 ERROR(sln, "%s:%s: error opening file -- %s\n", 307 caller, name, strerror(errno)); 308 res = RES_IO_ERR; 309 goto error; 310 } 311 } 312 313 exit: 314 return res; 315 error: 316 if(stream->intern_fp && stream->fp) CHK(fclose(stream->fp) == 0); 317 goto exit; 318 } 319 320 static res_T 321 create_tree 322 (struct sln_device* sln, 323 const char* caller, 324 struct sln_tree** out_tree) 325 { 326 struct sln_tree* tree = NULL; 327 res_T res = RES_OK; 328 ASSERT(sln && caller && out_tree); 329 330 tree = MEM_CALLOC(sln->allocator, 1, sizeof(struct sln_tree)); 331 if(!tree) { 332 ERROR(sln, "%s: could not allocate the tree data structure.\n", 333 caller); 334 res = RES_MEM_ERR; 335 goto error; 336 } 337 ref_init(&tree->ref); 338 SLN(device_ref_get(sln)); 339 tree->sln = sln; 340 darray_node_init(sln->allocator, &tree->nodes); 341 darray_vertex_init(sln->allocator, &tree->vertices); 342 343 exit: 344 *out_tree = tree; 345 return res; 346 error: 347 if(tree) { SLN(tree_ref_put(tree)); tree = NULL; } 348 goto exit; 349 } 350 351 static INLINE int 352 cmp_nu_vtx(const void* key, const void* item) 353 { 354 const float nu = *((const float*)key); 355 const struct sln_vertex* vtx = item; 356 if(nu < vtx->wavenumber) return -1; 357 if(nu > vtx->wavenumber) return +1; 358 return 0; 359 } 360 361 static void 362 release_tree(ref_T* ref) 363 { 364 struct sln_tree* tree = CONTAINER_OF(ref, struct sln_tree, ref); 365 struct sln_device* sln = NULL; 366 ASSERT(ref); 367 sln = tree->sln; 368 darray_node_release(&tree->nodes); 369 darray_vertex_release(&tree->vertices); 370 if(tree->args.lines) SHTR(line_list_ref_put(tree->args.lines)); 371 if(tree->args.metadata) SHTR(isotope_metadata_ref_put(tree->args.metadata)); 372 MEM_RM(sln->allocator, tree); 373 SLN(device_ref_put(sln)); 374 } 375 376 /******************************************************************************* 377 * Exported symbols 378 ******************************************************************************/ 379 res_T 380 sln_tree_create 381 (struct sln_device* device, 382 const struct sln_tree_create_args* args, 383 struct sln_tree** out_tree) 384 { 385 struct sln_tree* tree = NULL; 386 res_T res = RES_OK; 387 388 if(!device || !out_tree) { res = RES_BAD_ARG; goto error; } 389 res = check_sln_tree_create_args(device, FUNC_NAME, args); 390 if(res != RES_OK) goto error; 391 392 res = create_tree(device, FUNC_NAME, &tree); 393 if(res != RES_OK) goto error; 394 SHTR(line_list_ref_get(args->lines)); 395 SHTR(isotope_metadata_ref_get(args->metadata)); 396 tree->args = *args; 397 398 res = tree_build(tree); 399 if(res != RES_OK) goto error; 400 401 exit: 402 if(out_tree) *out_tree = tree; 403 return res; 404 error: 405 if(tree) { SLN(tree_ref_put(tree)); tree = NULL; } 406 goto exit; 407 } 408 409 res_T 410 sln_tree_read 411 (struct sln_device* sln, 412 const struct sln_tree_read_args* args, 413 struct sln_tree** out_tree) 414 { 415 struct shtr_line_list_read_args rlines_args = SHTR_LINE_LIST_READ_ARGS_NULL; 416 struct stream stream = STREAM_NULL; 417 struct sln_tree* tree = NULL; 418 size_t n = 0; 419 int version = 0; 420 res_T res = RES_OK; 421 422 if(!sln || !out_tree) { res = RES_BAD_ARG; goto error; } 423 res = check_sln_tree_read_args(sln, FUNC_NAME, args); 424 if(res != RES_OK) goto error; 425 426 res = create_tree(sln, FUNC_NAME, &tree); 427 if(res != RES_OK) goto error; 428 429 res = stream_init(sln, FUNC_NAME, args->filename, args->file, "r", &stream); 430 if(res != RES_OK) goto error; 431 432 #define READ(Var, Nb) { \ 433 if(fread((Var), sizeof(*(Var)), (Nb), stream.fp) != (Nb)) { \ 434 if(feof(stream.fp)) { \ 435 res = RES_BAD_ARG; \ 436 } else if(ferror(stream.fp)) { \ 437 res = RES_IO_ERR; \ 438 } else { \ 439 res = RES_UNKNOWN_ERR; \ 440 } \ 441 goto error; \ 442 } \ 443 } (void)0 444 READ(&version, 1); 445 if(version != SLN_TREE_VERSION) { 446 ERROR(sln, 447 "%s: unexpected tree version %d. Expecting a tree in version %d.\n", 448 FUNC_NAME, version, SLN_TREE_VERSION); 449 res = RES_BAD_ARG; 450 goto error; 451 } 452 453 READ(&n, 1); 454 if((res = darray_node_resize(&tree->nodes, n)) != RES_OK) goto error; 455 READ(darray_node_data_get(&tree->nodes), n); 456 457 READ(&n, 1); 458 if((res = darray_vertex_resize(&tree->vertices, n)) != RES_OK) goto error; 459 READ(darray_vertex_data_get(&tree->vertices), n); 460 461 READ(&tree->args.line_profile, 1); 462 READ(&tree->args.molecules, 1); 463 READ(&tree->args.pressure, 1); 464 READ(&tree->args.temperature, 1); 465 READ(&tree->args.nvertices_hint, 1); 466 READ(&tree->args.mesh_decimation_err, 1); 467 READ(&tree->args.mesh_type, 1); 468 #undef READ 469 470 res = shtr_isotope_metadata_create_from_stream 471 (args->shtr, stream.fp, &tree->args.metadata); 472 if(res != RES_OK) goto error; 473 474 rlines_args.file = stream.fp; 475 res = shtr_line_list_read(args->shtr, &rlines_args, &tree->args.lines); 476 if(res != RES_OK) goto error; 477 478 exit: 479 stream_release(&stream); 480 if(out_tree) *out_tree = tree; 481 return res; 482 error: 483 if(tree) { SLN(tree_ref_put(tree)); tree = NULL; } 484 if(sln) { 485 ERROR(sln, "%s: error loading the tree structure -- %s.\n", 486 FUNC_NAME, res_to_cstr(res)); 487 } 488 goto exit; 489 } 490 491 res_T 492 sln_tree_ref_get(struct sln_tree* tree) 493 { 494 if(!tree) return RES_BAD_ARG; 495 ref_get(&tree->ref); 496 return RES_OK; 497 } 498 499 res_T 500 sln_tree_ref_put(struct sln_tree* tree) 501 { 502 if(!tree) return RES_BAD_ARG; 503 ref_put(&tree->ref, release_tree); 504 return RES_OK; 505 } 506 507 res_T 508 sln_tree_get_desc(const struct sln_tree* tree, struct sln_tree_desc* desc) 509 { 510 if(!tree || !desc) return RES_BAD_ARG; 511 desc->max_nlines_per_leaf = 1; 512 desc->mesh_decimation_err = tree->args.mesh_decimation_err; 513 desc->mesh_type = tree->args.mesh_type; 514 desc->line_profile = tree->args.line_profile; 515 desc->nnodes = darray_node_size_get(&tree->nodes); 516 desc->nvertices = darray_vertex_size_get(&tree->vertices); 517 return RES_OK; 518 } 519 520 const struct sln_node* 521 sln_tree_get_root(const struct sln_tree* tree) 522 { 523 ASSERT(tree); 524 if(darray_node_size_get(&tree->nodes)) { 525 return darray_node_cdata_get(&tree->nodes); 526 } else { 527 return NULL; 528 } 529 } 530 531 int 532 sln_node_is_leaf(const struct sln_node* node) 533 { 534 ASSERT(node); 535 return node->offset == 0; 536 } 537 538 const struct sln_node* 539 sln_node_get_child(const struct sln_node* node, const unsigned ichild) 540 { 541 ASSERT(node && ichild <= 1); 542 return node + node->offset + ichild; 543 } 544 545 size_t 546 sln_node_get_lines_count(const struct sln_node* node) 547 { 548 ASSERT(node); 549 return node->range[1] - node->range[0] + 1/*Both boundaries are inclusives*/; 550 } 551 552 res_T 553 sln_node_get_line 554 (const struct sln_tree* tree, 555 const struct sln_node* node, 556 const size_t iline, 557 struct shtr_line* line) 558 { 559 if(!tree || !node || iline > sln_node_get_lines_count(node)) 560 return RES_BAD_ARG; 561 562 return shtr_line_list_at 563 (tree->args.lines, node->range[0] + iline, line); 564 } 565 566 res_T 567 sln_node_get_mesh 568 (const struct sln_tree* tree, 569 const struct sln_node* node, 570 struct sln_mesh* mesh) 571 { 572 if(!tree || !node || !mesh) return RES_BAD_ARG; 573 mesh->vertices = darray_vertex_cdata_get(&tree->vertices) + node->ivertex; 574 mesh->nvertices = node->nvertices; 575 return RES_OK; 576 } 577 578 double 579 sln_node_eval 580 (const struct sln_tree* tree, 581 const struct sln_node* node, 582 const double nu) 583 { 584 double ka = 0; 585 size_t iline; 586 ASSERT(tree && node); 587 588 FOR_EACH(iline, node->range[0], node->range[1]+1) { 589 struct line line = LINE_NULL; 590 res_T res = RES_OK; 591 592 res = line_setup(tree, iline, &line); 593 if(res != RES_OK) { 594 WARN(tree->sln, "%s: could not setup the line %lu-- %s\n", 595 FUNC_NAME, iline, res_to_cstr(res)); 596 } 597 598 ka += line_value(tree, &line, nu); 599 } 600 return ka; 601 } 602 603 double 604 sln_mesh_eval(const struct sln_mesh* mesh, const double wavenumber) 605 { 606 const struct sln_vertex* vtx0 = NULL; 607 const struct sln_vertex* vtx1 = NULL; 608 const float nu = (float)wavenumber; 609 size_t n; /* #vertices */ 610 float u; /* Linear interpolation paraeter */ 611 ASSERT(mesh && mesh->nvertices); 612 613 n = mesh->nvertices; 614 615 /* Handle special cases */ 616 if(n == 1) return mesh->vertices[0].ka; 617 if(nu <= mesh->vertices[0].wavenumber) return mesh->vertices[0].ka; 618 if(nu >= mesh->vertices[n-1].wavenumber) return mesh->vertices[n-1].ka; 619 620 /* Dichotomic search of the mesh vertex whose wavenumber is greater than or 621 * equal to the submitted wavenumber 'nu' */ 622 vtx1 = search_lower_bound(&nu, mesh->vertices, n, sizeof(*vtx1), cmp_nu_vtx); 623 vtx0 = vtx1 - 1; 624 ASSERT(vtx1); /* A vertex is necessary found ...*/ 625 ASSERT(vtx1 > mesh->vertices); /* ... and it cannot be the first one */ 626 ASSERT(vtx0->wavenumber < nu && nu <= vtx1->wavenumber); 627 628 /* Compute the linear interpolation parameter */ 629 u = (nu - vtx0->wavenumber) / (vtx1->wavenumber - vtx0->wavenumber); 630 u = CLAMP(u, 0, 1); /* Handle numerical imprecisions */ 631 632 if(u == 0) return vtx0->ka; 633 if(u == 1) return vtx1->ka; 634 return u*(vtx1->ka - vtx0->ka) + vtx0->ka; 635 } 636 637 res_T 638 sln_tree_write 639 (const struct sln_tree* tree, 640 const struct sln_tree_write_args* args) 641 { 642 struct shtr_line_list_write_args wlines_args = SHTR_LINE_LIST_WRITE_ARGS_NULL; 643 struct stream stream = STREAM_NULL; 644 size_t nnodes, nverts; 645 res_T res = RES_OK; 646 647 if(!tree) { res = RES_BAD_ARG; goto error; } 648 res = check_sln_tree_write_args(tree->sln, FUNC_NAME, args); 649 if(res != RES_OK) goto error; 650 651 res = stream_init 652 (tree->sln, FUNC_NAME, args->filename, args->file, "w", &stream); 653 if(res != RES_OK) goto error; 654 655 #define WRITE(Var, Nb) { \ 656 if(fwrite((Var), sizeof(*(Var)), (Nb), stream.fp) != (Nb)) { \ 657 ERROR(tree->sln, "%s:%s: error writing the tree -- %s\n", \ 658 FUNC_NAME, stream.name, strerror(errno)); \ 659 res = RES_IO_ERR; \ 660 goto error; \ 661 } \ 662 } (void)0 663 WRITE(&SLN_TREE_VERSION, 1); 664 665 nnodes = darray_node_size_get(&tree->nodes); 666 WRITE(&nnodes, 1); 667 WRITE(darray_node_cdata_get(&tree->nodes), nnodes); 668 669 nverts = darray_vertex_size_get(&tree->vertices); 670 WRITE(&nverts, 1); 671 WRITE(darray_vertex_cdata_get(&tree->vertices), nverts); 672 673 WRITE(&tree->args.line_profile, 1); 674 WRITE(&tree->args.molecules, 1); 675 WRITE(&tree->args.pressure, 1); 676 WRITE(&tree->args.temperature, 1); 677 WRITE(&tree->args.nvertices_hint, 1); 678 WRITE(&tree->args.mesh_decimation_err, 1); 679 WRITE(&tree->args.mesh_type, 1); 680 #undef WRITE 681 682 res = shtr_isotope_metadata_write(tree->args.metadata, stream.fp); 683 if(res != RES_OK) goto error; 684 685 wlines_args.file = stream.fp; 686 res = shtr_line_list_write(tree->args.lines, &wlines_args); 687 if(res != RES_OK) goto error; 688 689 exit: 690 stream_release(&stream); 691 return res; 692 error: 693 goto exit; 694 }