sln_slab.c (14943B)
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 /* Check that the probability is valid */ 315 static INLINE res_T 316 check_proba(const struct cmd* cmd, const double proba) 317 { 318 if(0 <= proba && proba <= 1) return RES_OK; 319 320 if(cmd->args.verbose >= 1) { 321 fprintf(stderr, "error: invalid probability %g\n", proba); 322 } 323 return RES_BAD_ARG; 324 } 325 326 static FINLINE double /* [W/m^2/sr/cm^-1] */ 327 planck 328 (const double nu/* [cm^-1] */, 329 const double T/* [K] */) 330 { 331 const double lambda = WAVENUMBER_TO_WAVELENGTH(nu); /* [m] */ 332 const double planck1 = sbb_planck_monochromatic(lambda, T); /* [W/m^2/sr/m] */ 333 const double planck2 = planck1 * (1.0e-2/(nu*nu)); /* [W/m^2/sr/cm^-1] */ 334 ASSERT(T >= 0); 335 return planck2; 336 } 337 338 static INLINE const char* 339 estimate_cstr(const enum estimate estimate) 340 { 341 const char* cstr = NULL; 342 switch(estimate) { 343 case TRANSMISSIVITY: cstr="transmissivity"; break; 344 case EMISSIVITY: cstr="emissivity"; break; 345 default: FATAL("Unreachable code\n"); break; 346 } 347 return cstr; 348 } 349 350 static res_T 351 realisation 352 (const struct cmd* cmd, 353 struct ssp_rng* rng, 354 double out_weights[ESTIMATE_COUNT__]) 355 { 356 /* Acceleration structure */ 357 struct sln_tree_desc tree_desc = SLN_TREE_DESC_NULL; 358 const struct sln_node* root = NULL; 359 struct sln_mesh mesh = SLN_MESH_NULL; 360 361 /* Null collisions */ 362 double ka_max = 0; 363 double dst_remain = 0; 364 size_t ncollisions = 0; /* Number of null collisions */ 365 366 /* Miscellaneous */ 367 double w[ESTIMATE_COUNT__] = {0, 0}; /* Monte Carlo weight */ 368 double nu = 0; /* Sampled wavenumber [cm^-1] */ 369 double nu_range = 0; /* Spectral range [cm^-1] */ 370 int i = 0; 371 res_T res = RES_OK; 372 373 ASSERT(cmd && rng && out_weights); /* Check pre-conditions */ 374 375 /* Precompute the spectral range */ 376 nu_range = cmd->args.spectral_range[1] - cmd->args.spectral_range[0]; 377 378 /* Initialize the total distance to traverse with the thickness of the slab */ 379 dst_remain = cmd->args.thickness; 380 381 /* Uniformly sample the spectral dimension */ 382 nu = ssp_rng_uniform_double 383 (rng, cmd->args.spectral_range[0], cmd->args.spectral_range[1]); 384 385 SLN(tree_get_desc(cmd->tree, &tree_desc)); 386 387 /* Retrieve the ka_max of the spectrum at the sampled nu */ 388 root = sln_tree_get_root(cmd->tree); 389 SLN(node_get_mesh(cmd->tree, root, &mesh)); 390 ka_max = sln_mesh_eval(&mesh, nu); 391 392 for(ncollisions=0; ; ++ncollisions) { 393 const struct sln_node* leaf = NULL; 394 double dst = 0; /* Sampled distance */ 395 double proba_abs = 0; /* Probability of absorption */ 396 double leaf_proba = 0; /* Probability of sampling the line */ 397 double leaf_ka = 0; /* Value of the line */ 398 double r = 0; /* Random number */ 399 400 dst = ssp_ran_exp(rng, ka_max); /* Sample a traversal distance */ 401 402 if(dst > dst_remain) { /* No absorption in the slab */ 403 w[TRANSMISSIVITY] = 1.0; 404 w[EMISSIVITY] = 0.0; 405 break; 406 } 407 408 /* Importance sampling of a line */ 409 leaf = sln_node_sample_leaf(cmd->tree, root, nu, rng, &leaf_proba); 410 if(!leaf) { res = RES_BAD_ARG; goto error; } 411 412 /* Evaluate the value of the line and compute the probability of being 413 * absorbed by it */ 414 leaf_ka = sln_node_eval(cmd->tree, leaf, nu); 415 proba_abs = leaf_ka / (leaf_proba*ka_max); 416 if((res = check_proba(cmd, proba_abs)) != RES_OK) goto error; 417 418 r = ssp_rng_canonical(rng); 419 if(r < proba_abs) { /* An absorption occurs */ 420 w[TRANSMISSIVITY] = 0.0; 421 w[EMISSIVITY] = planck(nu, tree_desc.temperature)*nu_range; /*[W/m^2/sr]*/ 422 break; 423 } 424 425 dst_remain -= dst; /* This was a null transition. Go on */ 426 } 427 428 exit: 429 FOR_EACH(i, 0, ESTIMATE_COUNT__) out_weights[i] = w[i]; 430 return res; 431 error: 432 FOR_EACH(i, 0, ESTIMATE_COUNT__) w[i] = NaN; 433 goto exit; 434 } 435 436 static res_T 437 cmd_run(const struct cmd* cmd) 438 { 439 /* Random Number Generator */ 440 struct ssp_rng** rngs = NULL; 441 442 /* Monte Carlo */ 443 struct accum accum[ESTIMATE_COUNT__] = {0}; 444 int64_t i = 0; /* Index of the realisation */ 445 size_t nrejects = 0; /* Number of rejected realisations */ 446 447 /* Progress */ 448 size_t nrealisations = 0; 449 size_t realisation_done = 0; 450 int progress = 0; 451 int progress_pcent = 10; 452 453 res_T res = RES_OK; 454 ASSERT(cmd); 455 456 res = create_per_thread_rngs(cmd, &rngs); 457 if(res != RES_OK) goto error; 458 459 #define PROGRESS_MSG "Solving: %3d%%\n" 460 if(cmd->args.verbose >= 3) fprintf(stderr, PROGRESS_MSG, progress); 461 462 nrealisations = cmd->args.nrealisations; 463 464 omp_set_num_threads((int)cmd->nthreads); 465 466 #pragma omp parallel for schedule(static) 467 for(i = 0; i < (int64_t)nrealisations; ++i) { 468 double w[ESTIMATE_COUNT__] = {0}; /* Monte Carlo weights */ 469 const int ithread = omp_get_thread_num(); 470 int pcent = 0; 471 res_T res_realisation = RES_OK; 472 473 res_realisation = realisation(cmd, rngs[ithread], w); 474 475 #pragma omp critical 476 { 477 /* Update the Monte Carlo accumulator */ 478 if(res_realisation == RES_OK) { 479 int iestim = 0; 480 FOR_EACH(iestim, 0, ESTIMATE_COUNT__) { 481 accum[iestim].sum += w[iestim]; 482 accum[iestim].sum2 += w[iestim]*w[iestim]; 483 accum[iestim].count += 1; 484 } 485 } 486 487 if(cmd->args.verbose >= 3) { 488 /* Update progress */ 489 realisation_done += 1; 490 pcent = (int)((double)realisation_done*100.0/(double)nrealisations+0.5); 491 if(pcent/progress_pcent > progress/progress_pcent) { 492 progress = pcent; 493 fprintf(stderr, PROGRESS_MSG, progress); 494 } 495 } 496 } 497 } 498 499 #undef PROGRESS_MSG 500 501 nrejects = nrealisations - accum[0].count; 502 503 FOR_EACH(i, 0, ESTIMATE_COUNT__) { 504 const double E = accum[i].sum / (double)accum[i].count; 505 const double V = accum[i].sum2 / (double)accum[i].count - E*E; 506 const double SE = sqrt(V/(double)accum[i].count); 507 508 /* Assume that the number of realisations is the same for all estimates */ 509 ASSERT(accum[i].count == accum[0].count); 510 511 printf("%-16s: %e +/- %e; %lu\n", estimate_cstr(i), E, SE, (unsigned long)nrejects); 512 } 513 514 exit: 515 delete_per_thread_rngs(cmd, rngs); 516 return res; 517 error: 518 goto exit; 519 } 520 521 /******************************************************************************* 522 * Main function 523 ******************************************************************************/ 524 int 525 main(int argc, char** argv) 526 { 527 struct args args = ARGS_DEFAULT; 528 struct cmd cmd = CMD_NULL; 529 int err = 0; 530 res_T res = RES_OK; 531 532 if((res = args_init(&args, argc, argv)) != RES_OK) goto error; 533 if(args.quit) goto exit; 534 535 if((res = cmd_init(&cmd, &args)) != RES_OK) goto error; 536 if((res = cmd_run(&cmd)) != RES_OK) goto error; 537 538 exit: 539 cmd_release(&cmd); 540 CHK(mem_allocated_size() == 0); 541 return err; 542 error: 543 err = 1; 544 goto exit; 545 }