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