sln_slab.c (10454B)
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 fprintf(stderr, 228 "Error creating the list of per thread RNG -- %s\n", 229 res_to_cstr(res)); 230 if(rngs) delete_per_thread_rngs(cmd, rngs); 231 rngs = NULL; 232 goto exit; 233 } 234 235 static void 236 cmd_release(struct cmd* cmd) 237 { 238 ASSERT(cmd); 239 if(cmd->tree) SLN(tree_ref_put(cmd->tree)); 240 } 241 242 static res_T 243 cmd_init(struct cmd* cmd, const struct args* args) 244 { 245 /* Star Line */ 246 struct sln_device_create_args sln_args = SLN_DEVICE_CREATE_ARGS_DEFAULT; 247 struct sln_tree_read_args tree_args = SLN_TREE_READ_ARGS_NULL; 248 struct sln_device* sln = NULL; 249 250 /* Star HITRAN */ 251 struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT; 252 struct shtr* shtr = NULL; 253 struct shtr_isotope_metadata* molparams = NULL; 254 struct shtr_line_list* lines = NULL; 255 256 /* Miscellaneous */ 257 unsigned nthreads_max = 0; 258 res_T res = RES_OK; 259 260 ASSERT(cmd && args); 261 262 *cmd = CMD_NULL; 263 264 shtr_args.verbose = args->verbose; 265 res = shtr_create(&shtr_args, &shtr); 266 if(res != RES_OK) goto error; 267 268 res = shtr_isotope_metadata_load(shtr, args->molparams, &molparams); 269 if(res != RES_OK) goto error; 270 271 res = load_lines(shtr, args, &lines); 272 if(res != RES_OK) goto error; 273 274 sln_args.verbose = args->verbose; 275 res = sln_device_create(&sln_args, &sln); 276 if(res != RES_OK) goto error; 277 278 tree_args.metadata = molparams; 279 tree_args.lines = lines; 280 tree_args.filename = args->tree; 281 res = sln_tree_read(sln, &tree_args, &cmd->tree); 282 if(res != RES_OK) goto error; 283 284 nthreads_max = (unsigned)MMAX(omp_get_max_threads(), omp_get_num_procs()); 285 cmd->args = *args; 286 cmd->nthreads = MMIN(cmd->args.nthreads_hint, nthreads_max); 287 288 exit: 289 if(sln) SLN(device_ref_put(sln)); 290 if(shtr) SHTR(ref_put(shtr)); 291 if(molparams) SHTR(isotope_metadata_ref_put(molparams)); 292 if(lines) SHTR(line_list_ref_put(lines)); 293 return res; 294 error: 295 cmd_release(cmd); 296 *cmd = CMD_NULL; 297 goto exit; 298 } 299 300 static double 301 realisation(const struct cmd* cmd, struct ssp_rng* rng) 302 { 303 ASSERT(cmd && rng); 304 (void)cmd; /* Avoid unused variable warning */ 305 306 if(ssp_rng_canonical(rng) < 0.5) { 307 return 2.0; 308 } else { 309 return 8.0; 310 } 311 } 312 313 static res_T 314 cmd_run(const struct cmd* cmd) 315 { 316 /* Random Number Generator */ 317 struct ssp_rng** rngs = NULL; 318 319 /* Monte Carlo */ 320 struct accum accum = ACCUM_NULL; 321 double E = 0; /* Expected value */ 322 double V = 0; /* Variance */ 323 double SE = 0; /* Standard Error */ 324 int64_t i = 0; /* Index of the realisation */ 325 326 /* Progress */ 327 size_t nrealisations = 0; 328 size_t realisation_done = 0; 329 int progress = 0; 330 int progress_pcent = 10; 331 332 res_T res = RES_OK; 333 ASSERT(cmd); 334 335 res = create_per_thread_rngs(cmd, &rngs); 336 if(res != RES_OK) goto error; 337 338 if(cmd->args.verbose) fprintf(stderr, "%3d%%\n", progress); 339 340 nrealisations = cmd->args.nrealisations; 341 342 #pragma omp parallel for schedule(static) 343 for(i = 0; i < (int64_t)nrealisations; ++i) { 344 const int ithread = omp_get_thread_num(); 345 int pcent = 0; 346 double w = 0; 347 348 w = realisation(cmd, rngs[ithread]); 349 350 #pragma omp critical 351 { 352 /* Update the Monte Carlo accumulator */ 353 accum.sum += w; 354 accum.sum2 += w*w; 355 accum.count += 1; 356 357 if(cmd->args.verbose) { 358 /* Update progress */ 359 realisation_done += 1; 360 pcent = (int)((double)realisation_done*100.0/(double)nrealisations+0.5); 361 if(pcent/progress_pcent > progress/progress_pcent) { 362 progress = pcent; 363 fprintf(stderr, "%3d%%\n", progress); 364 } 365 } 366 } 367 } 368 369 E = accum.sum / (double)accum.count; 370 V = accum.sum2 / (double)accum.count - E*E; 371 SE = sqrt(V/(double)accum.count); 372 373 printf("%g %g\n", E, SE); 374 375 exit: 376 delete_per_thread_rngs(cmd, rngs); 377 return res; 378 error: 379 goto exit; 380 } 381 382 /******************************************************************************* 383 * Main function 384 ******************************************************************************/ 385 int 386 main(int argc, char** argv) 387 { 388 struct args args = ARGS_DEFAULT; 389 struct cmd cmd = CMD_NULL; 390 int err = 0; 391 res_T res = RES_OK; 392 393 if((res = args_init(&args, argc, argv)) != RES_OK) goto error; 394 if(args.quit) goto exit; 395 396 if((res = cmd_init(&cmd, &args)) != RES_OK) goto error; 397 if((res = cmd_run(&cmd)) != RES_OK) goto error; 398 399 exit: 400 cmd_release(&cmd); 401 CHK(mem_allocated_size() == 0); 402 return err; 403 error: 404 err = 1; 405 goto exit; 406 }