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