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