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