sln_slab.c (18135B)
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 #define _POSIX_C_SOURCE 200112L /* getopt */ 20 21 #include "sln.h" 22 23 #include <star/shtr.h> 24 #include <star/sbb.h> 25 #include <star/ssp.h> 26 27 #include <rsys/cstr.h> 28 #include <rsys/mem_allocator.h> 29 #include <rsys/str.h> 30 31 #include <omp.h> 32 33 #include <unistd.h> /* getopt */ 34 35 enum estimate { 36 TRANSMISSIVITY, 37 EMISSIVITY, 38 ESTIMATE_COUNT__ 39 }; 40 41 #define WAVENUMBER_TO_WAVELENGTH(Nu/* [cm^-1] */) (1.e-2/(Nu))/*[m]*/ 42 43 struct args { 44 const char* tree; /* Acceleration structure */ 45 const char* molparams; 46 const char* lines; 47 48 double thickness; /* Thickness of the slab [m] */ 49 50 double spectral_range[2]; /* [cm^-1]^2 */ 51 52 unsigned long nrealisations; /* Number of Monte Carlo realisations */ 53 54 /* Miscellaneous */ 55 unsigned nthreads_hint; /* Hint on the number of threads to use */ 56 int lines_in_shtr_format; 57 int verbose; 58 int quit; 59 }; 60 #define ARGS_DEFAULT__ {NULL,NULL,NULL,1,{0,DBL_MAX},10000,UINT_MAX,0,0,0} 61 static const struct args ARGS_DEFAULT = ARGS_DEFAULT__; 62 63 struct cmd { 64 struct args args; 65 66 struct sln_tree* tree; 67 unsigned nthreads; 68 }; 69 #define CMD_NULL__ {0} 70 static const struct cmd CMD_NULL = CMD_NULL__; 71 72 struct accum { 73 double sum; 74 double sum2; 75 size_t count; 76 }; 77 #define ACCUM_NULL__ {0} 78 79 /******************************************************************************* 80 * Helper functions 81 ******************************************************************************/ 82 static void 83 usage(FILE* stream) 84 { 85 fprintf(stream, 86 "usage: sln-slab [-hsv] [-n nrealisations] [-T thickness] [-t threads]\n" 87 " -S nu_min,nu_max -a accel_struct -m molparams -l lines\n"); 88 } 89 90 static res_T 91 parse_spectral_range(const char* str, double spectral_range[2]) 92 { 93 size_t len = 0; 94 res_T res = RES_OK; 95 ASSERT(str && spectral_range); 96 97 res = cstr_to_list_double(str, ',', spectral_range, &len, 2); 98 if(res == RES_OK && len < 2) res = RES_BAD_ARG; 99 100 return res; 101 } 102 103 static res_T 104 args_init(struct args* args, int argc, char** argv) 105 { 106 int opt = 0; 107 res_T res = RES_OK; 108 109 ASSERT(args); 110 111 *args = ARGS_DEFAULT; 112 113 while((opt = getopt(argc, argv, "a:hl:m:n:S:sT:t:v")) != -1) { 114 switch(opt) { 115 case 'a': args->tree = optarg; break; 116 case 'h': 117 usage(stdout); 118 args->quit = 1; 119 goto exit; 120 case 'l': args->lines = optarg; break; 121 case 'm': args->molparams = optarg; break; 122 case 'n': res = cstr_to_ulong(optarg, &args->nrealisations); break; 123 case 'S': res = parse_spectral_range(optarg, args->spectral_range); break; 124 case 's': args->lines_in_shtr_format = 1; break; 125 case 'T': 126 res = cstr_to_double(optarg, &args->thickness); 127 if(res == RES_OK && args->thickness <= 0) res = RES_BAD_ARG; 128 break; 129 case 't': 130 res = cstr_to_uint(optarg, &args->nthreads_hint); 131 if(res == RES_OK && args->nthreads_hint == 0) res = RES_BAD_ARG; 132 break; 133 case 'v': args->verbose += (args->verbose < 3); break; 134 default: res = RES_BAD_ARG; break; 135 } 136 if(res != RES_OK) { 137 if(optarg) { 138 fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n", 139 argv[0], optarg, opt); 140 } 141 goto error; 142 } 143 } 144 145 #define MANDATORY(Cond, Name, Opt) { \ 146 if(!(Cond)) { \ 147 fprintf(stderr, "%s: %s missing -- option '-%c'\n", argv[0], (Name), (Opt)); \ 148 res = RES_BAD_ARG; \ 149 goto error; \ 150 } \ 151 } (void)0 152 MANDATORY(args->molparams, "molparams", 'm'); 153 MANDATORY(args->lines, "line list", 'l'); 154 MANDATORY(args->tree, "acceleration structure", 's'); 155 #undef MANDATORY 156 157 exit: 158 return res; 159 error: 160 usage(stderr); 161 goto exit; 162 } 163 164 static res_T 165 load_lines 166 (struct shtr* shtr, 167 const struct args* args, 168 struct shtr_line_list** out_lines) 169 { 170 struct shtr_line_list* lines = NULL; 171 res_T res = RES_OK; 172 ASSERT(shtr && args && out_lines); 173 174 if(args->lines_in_shtr_format) { 175 struct shtr_line_list_read_args read_args = SHTR_LINE_LIST_READ_ARGS_NULL; 176 177 /* Loads lines from data serialized by the Star-HITRAN library */ 178 read_args.filename = args->lines; 179 res = shtr_line_list_read(shtr, &read_args, &lines); 180 if(res != RES_OK) goto error; 181 182 } else { 183 struct shtr_line_list_load_args load_args = SHTR_LINE_LIST_LOAD_ARGS_NULL; 184 185 /* Loads lines from a file in HITRAN format */ 186 load_args.filename = args->lines; 187 res = shtr_line_list_load(shtr, &load_args, &lines); 188 if(res != RES_OK) goto error; 189 } 190 191 exit: 192 *out_lines = lines; 193 return res; 194 error: 195 if(lines) { SHTR(line_list_ref_put(lines)); lines = NULL; } 196 goto exit; 197 } 198 199 static void 200 delete_per_thread_rngs(const struct cmd* cmd, struct ssp_rng* rngs[]) 201 { 202 unsigned i = 0; 203 ASSERT(cmd && rngs); 204 205 FOR_EACH(i, 0, cmd->nthreads) { 206 if(rngs[i]) SSP(rng_ref_put(rngs[i])); 207 } 208 mem_rm(rngs); 209 } 210 211 static res_T 212 create_per_thread_rngs(const struct cmd* cmd, struct ssp_rng** out_rngs[]) 213 { 214 struct ssp_rng_proxy* proxy = NULL; 215 struct ssp_rng** rngs = NULL; 216 size_t i = 0; 217 res_T res = RES_OK; 218 ASSERT(cmd); 219 220 rngs = mem_calloc(cmd->nthreads, sizeof(*rngs)); 221 if(!rngs) { res = RES_MEM_ERR; goto error; } 222 223 res = ssp_rng_proxy_create(NULL, SSP_RNG_THREEFRY, cmd->nthreads, &proxy); 224 if(res != RES_OK) goto error; 225 226 FOR_EACH(i, 0, cmd->nthreads) { 227 res = ssp_rng_proxy_create_rng(proxy, i, &rngs[i]); 228 if(res != RES_OK) goto error; 229 } 230 231 exit: 232 *out_rngs = rngs; 233 if(proxy) SSP(rng_proxy_ref_put(proxy)); 234 return res; 235 error: 236 if(cmd->args.verbose >= 1) { 237 fprintf(stderr, 238 "Error creating the list of per thread RNG -- %s\n", 239 res_to_cstr(res)); 240 } 241 if(rngs) delete_per_thread_rngs(cmd, rngs); 242 rngs = NULL; 243 goto exit; 244 } 245 246 static void 247 cmd_release(struct cmd* cmd) 248 { 249 ASSERT(cmd); 250 if(cmd->tree) SLN(tree_ref_put(cmd->tree)); 251 } 252 253 static res_T 254 cmd_init(struct cmd* cmd, const struct args* args) 255 { 256 /* Star Line */ 257 struct sln_device_create_args sln_args = SLN_DEVICE_CREATE_ARGS_DEFAULT; 258 struct sln_tree_read_args tree_args = SLN_TREE_READ_ARGS_NULL; 259 struct sln_device* sln = NULL; 260 261 /* Star HITRAN */ 262 struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT; 263 struct shtr* shtr = NULL; 264 struct shtr_isotope_metadata* molparams = NULL; 265 struct shtr_line_list* lines = NULL; 266 267 /* Miscellaneous */ 268 unsigned nthreads_max = 0; 269 res_T res = RES_OK; 270 271 ASSERT(cmd && args); 272 273 *cmd = CMD_NULL; 274 275 shtr_args.verbose = args->verbose; 276 res = shtr_create(&shtr_args, &shtr); 277 if(res != RES_OK) goto error; 278 279 res = shtr_isotope_metadata_load(shtr, args->molparams, &molparams); 280 if(res != RES_OK) goto error; 281 282 res = load_lines(shtr, args, &lines); 283 if(res != RES_OK) goto error; 284 285 sln_args.verbose = args->verbose; 286 res = sln_device_create(&sln_args, &sln); 287 if(res != RES_OK) goto error; 288 289 tree_args.metadata = molparams; 290 tree_args.lines = lines; 291 tree_args.filename = args->tree; 292 res = sln_tree_read(sln, &tree_args, &cmd->tree); 293 if(res != RES_OK) goto error; 294 295 nthreads_max = (unsigned)MMAX(omp_get_max_threads(), omp_get_num_procs()); 296 cmd->args = *args; 297 cmd->nthreads = MMIN(cmd->args.nthreads_hint, nthreads_max); 298 299 exit: 300 if(sln) SLN(device_ref_put(sln)); 301 if(shtr) SHTR(ref_put(shtr)); 302 if(molparams) SHTR(isotope_metadata_ref_put(molparams)); 303 if(lines) SHTR(line_list_ref_put(lines)); 304 return res; 305 error: 306 cmd_release(cmd); 307 *cmd = CMD_NULL; 308 goto exit; 309 } 310 311 static res_T 312 build_node_cumulative 313 (const struct cmd* cmd, 314 const struct sln_node* node, 315 const double nu, /* [cm^-1] */ 316 const char* path, /* String representing the path to the node */ 317 double probabilities[SLN_TREE_ARITY_MAX], 318 double cumulative[SLN_TREE_ARITY_MAX]) 319 { 320 struct sln_mesh mesh = SLN_MESH_NULL; 321 double ka = 0; 322 double cumul = 0; 323 unsigned i=0, n=0; 324 res_T res = RES_OK; 325 ASSERT(cmd && node && path && cumulative); 326 327 n = sln_node_get_child_count(cmd->tree, node); 328 ASSERT(n <= SLN_TREE_ARITY_MAX); 329 330 FOR_EACH(i, 0, n) { 331 const struct sln_node* child = sln_node_get_child(cmd->tree, node, i); 332 333 SLN(node_get_mesh(cmd->tree, child, &mesh)); 334 335 ka = sln_mesh_eval(&mesh, nu); 336 337 cumul += ka; 338 cumulative[i] = cumul; 339 probabilities[i] = ka; 340 } 341 342 /* Check the criterion of transition importance sampling, i.e. the value of 343 * the parent node is greater than or equal to the sum of the values of its 344 * children */ 345 SLN(node_get_mesh(cmd->tree, node, &mesh)); 346 ka = sln_mesh_eval(&mesh, nu); 347 if(ka < cumul) { 348 if(cmd->args.verbose >= 1) { 349 fprintf(stderr, 350 "error: ka < ka_{0} + ka_{1} + ... ka_{N-1}; " 351 "%e < %e; nu=%-21.20g cm^-1; node path: %c%s\n", 352 ka, cumul, nu, path[0] != '\0' ? '-' : ' ', path); 353 } 354 res = RES_BAD_ARG; 355 goto error; 356 } 357 358 /* Complete the probability calculation and normalize the cumulative */ 359 ASSERT(cumul != 0); 360 FOR_EACH(i, 0, n) probabilities[i] /= cumul; 361 FOR_EACH(i, 0, n-1) cumulative[i] /= cumul; 362 cumulative[n-1] = 1; /* Handle numerical uncertainty */ 363 364 exit: 365 return res; 366 error: 367 goto exit; 368 } 369 370 static const struct sln_node* /* NULL <=> an error occurs */ 371 sample_line 372 (const struct cmd* cmd, 373 struct ssp_rng* rng, 374 const double nu/*[cm^-1]*/, 375 double* out_proba) 376 { 377 double cumulative[SLN_TREE_ARITY_MAX]; 378 double probabilities[SLN_TREE_ARITY_MAX]; 379 struct str path; 380 const struct sln_node* node = NULL; 381 double proba = 1; /* Proba to sample the line */ 382 size_t depth = 0; 383 res_T res = RES_OK; 384 ASSERT(cmd && rng); 385 386 str_init(NULL, &path); 387 node = sln_tree_get_root(cmd->tree); 388 389 for(depth=0; !sln_node_is_leaf(node); ++depth) { 390 double r = 0; 391 unsigned ichild = 0; 392 393 res = build_node_cumulative 394 (cmd, node, nu, str_cget(&path), probabilities, cumulative); 395 if(res != RES_OK) goto error; 396 397 /* Sample a child with respect to its importance */ 398 r = ssp_rng_canonical(rng); 399 FOR_EACH(ichild, 0, SLN_TREE_ARITY_MAX) { 400 if(r < cumulative[ichild]) { 401 proba *= probabilities[ichild]; 402 node = sln_node_get_child(cmd->tree, node, ichild); 403 break; 404 } 405 406 /* Update the path string */ 407 res = str_append_printf(&path, " -c%d", ichild); 408 if(res != RES_OK) goto error; 409 } 410 ASSERT(ichild < SLN_TREE_ARITY_MAX); /* A node should have been sampled */ 411 } 412 413 exit: 414 str_release(&path); 415 *out_proba = proba; 416 return node; 417 error: 418 proba = NaN; 419 node = NULL; 420 goto exit; 421 } 422 423 static INLINE res_T 424 check_sampled_node(const struct cmd* cmd, const struct sln_node* node) 425 { 426 ASSERT(cmd); 427 (void) cmd; 428 429 if(!node) return RES_BAD_ARG; 430 431 #ifndef NDEBUG 432 { 433 /* Verify that the sampled node corresponds to a single line */ 434 struct sln_node_desc desc = SLN_NODE_DESC_NULL; 435 SLN(node_get_desc(cmd->tree, node, &desc)); 436 ASSERT(desc.nlines == 1); 437 } 438 #endif 439 440 return RES_OK; 441 } 442 443 /* Check that the probability is valid */ 444 static INLINE res_T 445 check_proba(const struct cmd* cmd, const double proba) 446 { 447 if(0 <= proba && proba <= 1) return RES_OK; 448 449 if(cmd->args.verbose >= 1) { 450 fprintf(stderr, "error: invalid probability %g\n", proba); 451 } 452 return RES_BAD_ARG; 453 } 454 455 static FINLINE double /* [W/m^2/sr/cm^-1] */ 456 planck 457 (const double nu/* [cm^-1] */, 458 const double T/* [K] */) 459 { 460 const double lambda = WAVENUMBER_TO_WAVELENGTH(nu); /* [m] */ 461 const double planck1 = sbb_planck_monochromatic(lambda, T); /* [W/m^2/sr/m] */ 462 const double planck2 = planck1 * (1.0e-2/(nu*nu)); /* [W/m^2/sr/cm^-1] */ 463 ASSERT(T >= 0); 464 return planck2; 465 } 466 467 static INLINE const char* 468 estimate_cstr(const enum estimate estimate) 469 { 470 const char* cstr = NULL; 471 switch(estimate) { 472 case TRANSMISSIVITY: cstr="transmissivity"; break; 473 case EMISSIVITY: cstr="emissivity"; break; 474 default: FATAL("Unreachable code\n"); break; 475 } 476 return cstr; 477 } 478 479 static res_T 480 realisation 481 (const struct cmd* cmd, 482 struct ssp_rng* rng, 483 double out_weights[ESTIMATE_COUNT__]) 484 { 485 /* Acceleration structure */ 486 struct sln_tree_desc tree_desc = SLN_TREE_DESC_NULL; 487 const struct sln_node* node = NULL; 488 struct sln_mesh mesh = SLN_MESH_NULL; 489 490 /* Null collisions */ 491 double ka_max = 0; 492 double dst_remain = 0; 493 size_t ncollisions = 0; /* Number of null collisions */ 494 495 /* Miscellaneous */ 496 double w[ESTIMATE_COUNT__] = {0, 0}; /* Monte Carlo weight */ 497 double nu = 0; /* Sampled wavenumber [cm^-1] */ 498 double nu_range = 0; /* Spectral range [cm^-1] */ 499 int i = 0; 500 res_T res = RES_OK; 501 502 ASSERT(cmd && rng && out_weights); /* Check pre-conditions */ 503 504 /* Precompute the spectral range */ 505 nu_range = cmd->args.spectral_range[1] - cmd->args.spectral_range[0]; 506 507 /* Initialize the total distance to traverse with the thickness of the slab */ 508 dst_remain = cmd->args.thickness; 509 510 /* Uniformly sample the spectral dimension */ 511 nu = ssp_rng_uniform_double 512 (rng, cmd->args.spectral_range[0], cmd->args.spectral_range[1]); 513 514 SLN(tree_get_desc(cmd->tree, &tree_desc)); 515 516 /* Retrieve the ka_max of the spectrum at the sampled nu */ 517 node = sln_tree_get_root(cmd->tree); 518 SLN(node_get_mesh(cmd->tree, node, &mesh)); 519 ka_max = sln_mesh_eval(&mesh, nu); 520 521 for(ncollisions=0; ; ++ncollisions) { 522 double dst = 0; /* Sampled distance */ 523 double proba_abs = 0; /* Probability of absorption */ 524 double line_proba = 0; /* Probability of sampling the line */ 525 double line_ha = 0; /* Value of the line */ 526 double r = 0; /* Random number */ 527 528 dst = ssp_ran_exp(rng, ka_max); /* Sample a traversal distance */ 529 530 if(dst > dst_remain) { /* No absorption in the slab */ 531 w[TRANSMISSIVITY] = 1.0; 532 w[EMISSIVITY] = 0.0; 533 break; 534 } 535 536 /* Importance sampling of a line */ 537 node = sample_line(cmd, rng, nu, &line_proba); 538 if((res = check_sampled_node(cmd, node)) != RES_OK) goto error; 539 540 /* Evaluate the value of the line and compute the probability of being 541 * absorbed by it */ 542 line_ha = sln_node_eval(cmd->tree, node, nu); 543 proba_abs = line_ha / (line_proba*ka_max); 544 if((res = check_proba(cmd, proba_abs)) != RES_OK) goto error; 545 546 r = ssp_rng_canonical(rng); 547 if(r < proba_abs) { /* An absorption occurs */ 548 w[TRANSMISSIVITY] = 0.0; 549 w[EMISSIVITY] = planck(nu, tree_desc.temperature)*nu_range; /*[W/m^2/sr]*/ 550 break; 551 } 552 553 dst_remain -= dst; /* This was a null transition. Go on */ 554 } 555 556 exit: 557 FOR_EACH(i, 0, ESTIMATE_COUNT__) out_weights[i] = w[i]; 558 return res; 559 error: 560 FOR_EACH(i, 0, ESTIMATE_COUNT__) w[i] = NaN; 561 goto exit; 562 } 563 564 static res_T 565 cmd_run(const struct cmd* cmd) 566 { 567 /* Random Number Generator */ 568 struct ssp_rng** rngs = NULL; 569 570 /* Monte Carlo */ 571 struct accum accum[ESTIMATE_COUNT__] = {0}; 572 int64_t i = 0; /* Index of the realisation */ 573 size_t nrejects = 0; /* Number of rejected realisations */ 574 575 /* Progress */ 576 size_t nrealisations = 0; 577 size_t realisation_done = 0; 578 int progress = 0; 579 int progress_pcent = 10; 580 581 res_T res = RES_OK; 582 ASSERT(cmd); 583 584 res = create_per_thread_rngs(cmd, &rngs); 585 if(res != RES_OK) goto error; 586 587 #define PROGRESS_MSG "Solving: %3d%%\n" 588 if(cmd->args.verbose >= 3) fprintf(stderr, PROGRESS_MSG, progress); 589 590 nrealisations = cmd->args.nrealisations; 591 592 omp_set_num_threads((int)cmd->nthreads); 593 594 #pragma omp parallel for schedule(static) 595 for(i = 0; i < (int64_t)nrealisations; ++i) { 596 double w[ESTIMATE_COUNT__] = {0}; /* Monte Carlo weights */ 597 const int ithread = omp_get_thread_num(); 598 int pcent = 0; 599 res_T res_realisation = RES_OK; 600 601 res_realisation = realisation(cmd, rngs[ithread], w); 602 603 #pragma omp critical 604 { 605 /* Update the Monte Carlo accumulator */ 606 if(res_realisation == RES_OK) { 607 int iestim = 0; 608 FOR_EACH(iestim, 0, ESTIMATE_COUNT__) { 609 accum[iestim].sum += w[iestim]; 610 accum[iestim].sum2 += w[iestim]*w[iestim]; 611 accum[iestim].count += 1; 612 } 613 } 614 615 if(cmd->args.verbose >= 3) { 616 /* Update progress */ 617 realisation_done += 1; 618 pcent = (int)((double)realisation_done*100.0/(double)nrealisations+0.5); 619 if(pcent/progress_pcent > progress/progress_pcent) { 620 progress = pcent; 621 fprintf(stderr, PROGRESS_MSG, progress); 622 } 623 } 624 } 625 } 626 627 #undef PROGRESS_MSG 628 629 nrejects = nrealisations - accum[0].count; 630 631 FOR_EACH(i, 0, ESTIMATE_COUNT__) { 632 const double E = accum[i].sum / (double)accum[i].count; 633 const double V = accum[i].sum2 / (double)accum[i].count - E*E; 634 const double SE = sqrt(V/(double)accum[i].count); 635 636 /* Assume that the number of realisations is the same for all estimates */ 637 ASSERT(accum[i].count == accum[0].count); 638 639 printf("%-16s: %e +/- %e; %lu\n", estimate_cstr(i), E, SE, (unsigned long)nrejects); 640 } 641 642 exit: 643 delete_per_thread_rngs(cmd, rngs); 644 return res; 645 error: 646 goto exit; 647 } 648 649 /******************************************************************************* 650 * Main function 651 ******************************************************************************/ 652 int 653 main(int argc, char** argv) 654 { 655 struct args args = ARGS_DEFAULT; 656 struct cmd cmd = CMD_NULL; 657 int err = 0; 658 res_T res = RES_OK; 659 660 if((res = args_init(&args, argc, argv)) != RES_OK) goto error; 661 if(args.quit) goto exit; 662 663 if((res = cmd_init(&cmd, &args)) != RES_OK) goto error; 664 if((res = cmd_run(&cmd)) != RES_OK) goto error; 665 666 exit: 667 cmd_release(&cmd); 668 CHK(mem_allocated_size() == 0); 669 return err; 670 error: 671 err = 1; 672 goto exit; 673 }