star-line

Structure for accelerating line importance sampling
git clone git://git.meso-star.fr/star-line.git
Log | Files | Refs | README | LICENSE

sln_build.c (14578B)


      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 200809L /* strtok_r */
     20 
     21 #include "sln.h"
     22 
     23 #include <rsys/cstr.h>
     24 #include <rsys/mem_allocator.h>
     25 #include <rsys/rsys.h>
     26 
     27 #include <stdio.h>
     28 #include <string.h>
     29 #include <unistd.h> /* getopt */
     30 
     31 enum line_list_format {
     32   LINE_LIST_HITRAN,
     33   LINE_LIST_SHTR,
     34   LINE_LIST_FORMAT_COUNT__
     35 };
     36 
     37 struct args {
     38   /* Spectroscopic parameters */
     39   const char* lines;
     40   const char* molparams;
     41   enum sln_line_profile line_profile;
     42 
     43   const char* output;
     44 
     45   /* Thermodynamic properties */
     46   const char* mixture;
     47   double pressure; /* [atm] */
     48   double temperature; /* [K] */
     49 
     50   /* Polyline */
     51   double mesh_decimation_err;
     52   unsigned nvertices_hint;
     53   enum sln_mesh_type mesh_type;
     54   int collapse_polylines;
     55 
     56   /* Miscellaneous */
     57   int quit;
     58   int verbose;
     59   unsigned arity; /* tree arity */
     60   unsigned leaf_nlines; /* Maximum number of lines per leaf */
     61   unsigned nthreads_hint; /* Advice on the number of threads to use */
     62   enum line_list_format line_format;
     63 };
     64 #define ARGS_DEFAULT__ {                                                       \
     65   /* Spectroscopic parameters */                                               \
     66   NULL, /* line list */                                                        \
     67   NULL, /* Isotopologue metadata */                                            \
     68   SLN_LINE_PROFILE_VOIGT, /* Line profile */                                   \
     69                                                                                \
     70   NULL, /* Output */                                                           \
     71                                                                                \
     72   /* Thermodynamic properties */                                               \
     73   NULL,                                                                        \
     74   -1, /* Pressure [atm] */                                                     \
     75   -1, /* Temperature [K] */                                                    \
     76                                                                                \
     77   /* Polyline */                                                               \
     78   0.01,                                                                        \
     79   16,                                                                          \
     80   SLN_MESH_UPPER,                                                              \
     81   0,                                                                           \
     82                                                                                \
     83   /* Miscellaneous */                                                          \
     84   0, /* Quit */                                                                \
     85   0, /* Verbose */                                                             \
     86   2, /* Tree arity */                                                          \
     87   1, /* Number of lines per leaf */                                            \
     88   (unsigned)(-1), /* #threads hint */                                          \
     89   LINE_LIST_HITRAN /* lines_format */                                          \
     90 }
     91 static const struct args ARGS_DEFAULT = ARGS_DEFAULT__;
     92 
     93 struct cmd {
     94   struct args args;
     95 
     96   struct sln_device* sln;
     97   struct sln_mixture* mixture;
     98 
     99   struct shtr* shtr;
    100   struct shtr_line_list* lines;
    101   struct shtr_isotope_metadata* molparam;
    102 
    103   FILE* output;
    104 };
    105 static const struct cmd CMD_NULL = {0};
    106 
    107 /*******************************************************************************
    108  * Helper functions
    109  ******************************************************************************/
    110 static void
    111 usage(FILE* stream)
    112 {
    113   fprintf(stream,
    114 "usage: sln-build [-chsv] [-a arity] [-e polyline_opt[:polyline_opt ...]]\n"
    115 "                 [-L leaf_nlines] [-l line_profile] [-o accel_struct]\n"
    116 "                 [-t thread_count] -P pressure -T temperature -m molparms\n"
    117 "                 -x mixture [lines]\n");
    118 }
    119 
    120 static res_T
    121 parse_line_profile(const char* str, enum sln_line_profile* profile)
    122 {
    123   res_T res = RES_OK;
    124   ASSERT(str && profile);
    125 
    126   if(!strcmp(str, "voigt")) {
    127     *profile = SLN_LINE_PROFILE_VOIGT;
    128 
    129   } else {
    130     fprintf(stderr, "invalid line profile '%s'\n", str);
    131     res = RES_BAD_ARG;
    132     goto error;
    133   }
    134 
    135 exit:
    136   return res;
    137 error:
    138   goto exit;
    139 }
    140 
    141 static res_T
    142 parse_mesh_type(const char* str, enum sln_mesh_type* type)
    143 {
    144   res_T res = RES_OK;
    145   ASSERT(str && type);
    146 
    147   if(!strcmp(str, "fit")) {
    148     *type = SLN_MESH_FIT;
    149   } else if(!strcmp(str, "upper")) {
    150     *type = SLN_MESH_UPPER;
    151   } else {
    152     fprintf(stderr, "invalid mesh type `%s'\n", str);
    153     res = RES_BAD_ARG;
    154     goto error;
    155   }
    156 
    157 exit:
    158   return res;
    159 error:
    160   goto exit;
    161 }
    162 
    163 static res_T
    164 parse_polyline_opt(const char* str, void* ptr)
    165 {
    166   enum { ERR, MESH, VCOUNT } opt;
    167   char buf[BUFSIZ];
    168 
    169   struct args* args = ptr;
    170 
    171   char* key = NULL;
    172   char* val = NULL;
    173   char* tk_ctx = NULL;
    174   res_T res = RES_OK;
    175 
    176   ASSERT(str && ptr);
    177 
    178   if(strlen(str)+1/*NULL char*/ > sizeof(buf)) {
    179     fprintf(stderr, "could not duplicate polyline option `%s'\n", str);
    180     res = RES_MEM_ERR;
    181     goto error;
    182   }
    183 
    184   strncpy(buf, str, sizeof(buf));
    185 
    186   key = strtok_r(buf, "=", &tk_ctx);
    187   val = strtok_r(NULL, "", &tk_ctx);
    188 
    189        if(!strcmp(key, "err")) opt = ERR;
    190   else if(!strcmp(key, "mesh")) opt = MESH;
    191   else if(!strcmp(key, "vcount")) opt = VCOUNT;
    192   else {
    193     fprintf(stderr, "invalid polyline option `%s'\n", key);
    194     res = RES_BAD_ARG;
    195     goto error;
    196   }
    197 
    198   switch(opt) {
    199     case ERR:
    200       res = cstr_to_double(val, &args->mesh_decimation_err);
    201       if(res == RES_OK && args->mesh_decimation_err < 0) res = RES_BAD_ARG;
    202       break;
    203     case MESH:
    204       res = parse_mesh_type(val, &args->mesh_type);
    205       break;
    206     case VCOUNT:
    207       res = cstr_to_uint(val, &args->nvertices_hint);
    208       break;
    209     default: FATAL("Unreachable code\n"); break;
    210   }
    211 
    212   if(res != RES_OK) {
    213     fprintf(stderr,
    214       "error while parsing the polyline option `%s' -- %s\n",
    215       str, res_to_cstr(res));
    216     goto error;
    217   }
    218 
    219 exit:
    220   return res;
    221 error:
    222   goto exit;
    223 }
    224 
    225 static res_T
    226 args_init(struct args* args, int argc, char** argv)
    227 {
    228   int opt = 0;
    229   res_T res = RES_OK;
    230 
    231   ASSERT(args);
    232 
    233   *args = ARGS_DEFAULT;
    234 
    235   while((opt = getopt(argc, argv, "a:ce:hL:l:m:o:P:sT:t:vx:")) != -1) {
    236     switch(opt) {
    237       case 'a':
    238         res = cstr_to_uint(optarg, &args->arity);
    239         if(res == RES_OK && args->arity < 2) res = RES_BAD_ARG;
    240         break;
    241       case 'c': args->collapse_polylines = 1; break;
    242       case 'e':
    243         res = cstr_parse_list(optarg, ':', parse_polyline_opt, args);
    244         break;
    245       case 'h':
    246         usage(stdout);
    247         args->quit = 1;
    248         goto exit;
    249       case 'L':
    250         res = cstr_to_uint(optarg, &args->leaf_nlines);
    251         if(res == RES_OK && args->leaf_nlines < 1) res = RES_BAD_ARG;
    252         break;
    253       case 'l':
    254         res = parse_line_profile(optarg, &args->line_profile);
    255         break;
    256       case 'o': args->output = optarg; break;
    257       case 'P':
    258         res = cstr_to_double(optarg, &args->pressure);
    259         if(res == RES_OK && args->pressure < 0) res = RES_BAD_ARG;
    260         break;
    261       case 's': args->line_format = LINE_LIST_SHTR; break;
    262       case 'T':
    263         res = cstr_to_double(optarg, &args->temperature);
    264         if(res == RES_OK && args->temperature < 0) res = RES_BAD_ARG;
    265         break;
    266       case 't':
    267         res = cstr_to_uint(optarg, &args->nthreads_hint);
    268         if(res == RES_OK && args->nthreads_hint < 1) res = RES_BAD_ARG;
    269         break;
    270       case 'm': args->molparams = optarg; break;
    271       case 'v': args->verbose += (args->verbose < 3); break;
    272       case 'x': args->mixture = optarg; break;
    273       default: res = RES_BAD_ARG; break;
    274     }
    275 
    276     if(res != RES_OK) {
    277       if(optarg) {
    278         fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n",
    279           argv[0], optarg, opt);
    280       }
    281       goto error;
    282     }
    283   }
    284 
    285   if(optind < argc) args->lines = argv[optind];
    286 
    287   #define MANDATORY(Cond, Name, Opt) { \
    288     if(!(Cond)) { \
    289       fprintf(stderr, "%s: %s missing -- option '-%c'\n", argv[0], (Name), (Opt)); \
    290       res = RES_BAD_ARG; \
    291       goto error; \
    292     } \
    293   } (void)0
    294   MANDATORY(args->pressure >= 0, "pressure", 'P');
    295   MANDATORY(args->temperature >= 0, "temperature", 'T');
    296   MANDATORY(args->molparams, "molparams", 'm');
    297   MANDATORY(args->mixture, "mixture", 'x');
    298   #undef MANDATORY
    299 
    300 exit:
    301   return res;
    302 error:
    303   usage(stderr);
    304   goto exit;
    305 }
    306 
    307 static res_T
    308 load_lines_hitran(struct cmd* cmd, const struct args* args)
    309 {
    310   struct shtr_line_list_load_args load_args = SHTR_LINE_LIST_LOAD_ARGS_NULL;
    311   ASSERT(cmd && args);
    312 
    313   if(args->lines != NULL) {
    314     load_args.filename = args->lines;
    315   } else {
    316     load_args.filename = "stdin";
    317     load_args.file = stdin;
    318   }
    319 
    320   return shtr_line_list_load(cmd->shtr, &load_args, &cmd->lines);
    321 }
    322 
    323 static res_T
    324 load_lines_shtr(struct cmd* cmd, const struct args* args)
    325 {
    326   struct shtr_line_list_read_args rlines_args = SHTR_LINE_LIST_READ_ARGS_NULL;
    327 
    328   ASSERT(cmd && args);
    329 
    330   rlines_args.filename = args->lines;
    331   return shtr_line_list_read(cmd->shtr, &rlines_args, &cmd->lines);
    332 }
    333 
    334 static res_T
    335 load_lines(struct cmd* cmd, const struct args* args)
    336 {
    337   res_T res = RES_OK;
    338   switch(args->line_format) {
    339     case LINE_LIST_HITRAN: res = load_lines_hitran(cmd, args); break;
    340     case LINE_LIST_SHTR: res = load_lines_shtr(cmd, args); break;
    341     default: FATAL("Unreachable code\n"); break;
    342   }
    343   return res;
    344 }
    345 
    346 static res_T
    347 setup_output(struct cmd* cmd, const struct args* args)
    348 {
    349   res_T res = RES_OK;
    350   ASSERT(cmd && args);
    351 
    352   if(!args->output) {
    353     cmd->output = stdout;
    354 
    355   } else {
    356     cmd->output = fopen(args->output, "w");
    357     if(!cmd->output) {
    358       fprintf(stderr, "error opening file `%s' -- %s\n",
    359         args->output, strerror(errno));
    360       res = RES_IO_ERR;
    361       goto error;
    362     }
    363   }
    364 
    365 exit:
    366   return res;
    367 error:
    368   if(cmd->output && cmd->output != stdout) CHK(fclose(cmd->output) == 0);
    369   goto exit;
    370 }
    371 
    372 static res_T
    373 setup_tree_mixture
    374   (const struct cmd* cmd,
    375    struct sln_molecule molecules[SHTR_MAX_MOLECULE_COUNT])
    376 {
    377   int i = 0;
    378   int n = 0;
    379   res_T res = RES_OK;
    380   ASSERT(cmd && molecules);
    381 
    382   n = sln_mixture_get_molecule_count(cmd->mixture);
    383   FOR_EACH(i, 0, n) {
    384     enum shtr_molecule_id id = sln_mixture_get_molecule_id(cmd->mixture, i);
    385     res = sln_mixture_get_molecule(cmd->mixture, i, molecules+id);
    386     if(res != RES_OK) goto error;
    387   }
    388 
    389 exit:
    390   return res;
    391 error:
    392   goto exit;
    393 }
    394 
    395 static void
    396 cmd_release(struct cmd* cmd)
    397 {
    398   ASSERT(cmd);
    399   if(cmd->sln) SLN(device_ref_put(cmd->sln));
    400   if(cmd->mixture) SLN(mixture_ref_put(cmd->mixture));
    401   if(cmd->shtr) SHTR(ref_put(cmd->shtr));
    402   if(cmd->lines) SHTR(line_list_ref_put(cmd->lines));
    403   if(cmd->molparam) SHTR(isotope_metadata_ref_put(cmd->molparam));
    404 }
    405 
    406 static res_T
    407 cmd_init(struct cmd* cmd, const struct args* args)
    408 {
    409   struct sln_device_create_args sln_args = SLN_DEVICE_CREATE_ARGS_DEFAULT;
    410   struct sln_mixture_load_args mixture_args = SLN_MIXTURE_LOAD_ARGS_NULL;
    411   struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT;
    412   res_T res = RES_OK;
    413 
    414   ASSERT(cmd && args);
    415 
    416   *cmd = CMD_NULL;
    417 
    418   shtr_args.verbose = args->verbose;
    419   res = shtr_create(&shtr_args, &cmd->shtr);
    420   if(res != RES_OK) goto error;
    421 
    422   sln_args.verbose = args->verbose;
    423   res = sln_device_create(&sln_args, &cmd->sln);
    424   if(res != RES_OK) goto error;
    425 
    426   res = shtr_isotope_metadata_load(cmd->shtr, args->molparams, &cmd->molparam);
    427   if(res != RES_OK) goto error;
    428 
    429   mixture_args.filename = args->mixture;
    430   mixture_args.molparam = cmd->molparam;
    431   res = sln_mixture_load(cmd->sln, &mixture_args, &cmd->mixture);
    432   if(res != RES_OK) goto error;
    433 
    434   res = setup_output(cmd, args);
    435   if(res != RES_OK) goto error;
    436 
    437   res = load_lines(cmd, args);
    438   if(res != RES_OK) goto error;
    439 
    440   cmd->args = *args;
    441 
    442 exit:
    443   return res;
    444 error:
    445   cmd_release(cmd);
    446   *cmd = CMD_NULL;
    447   goto exit;
    448 }
    449 
    450 static res_T
    451 cmd_run(struct cmd* cmd)
    452 {
    453   struct sln_tree_create_args tree_args = SLN_TREE_CREATE_ARGS_DEFAULT;
    454   struct sln_tree_write_args write_args = SLN_TREE_WRITE_ARGS_NULL;
    455   struct sln_tree* tree = NULL;
    456   res_T res = RES_OK;
    457 
    458   ASSERT(cmd);
    459 
    460   tree_args.metadata = cmd->molparam;
    461   tree_args.lines = cmd->lines;
    462   tree_args.pressure = cmd->args.pressure;
    463   tree_args.temperature = cmd->args.temperature;
    464   tree_args.nvertices_hint = cmd->args.nvertices_hint;
    465   tree_args.mesh_decimation_err = cmd->args.mesh_decimation_err;
    466   tree_args.mesh_type = cmd->args.mesh_type;
    467   tree_args.arity = cmd->args.arity;
    468   tree_args.leaf_nlines = cmd->args.leaf_nlines;
    469   tree_args.collapse_polylines = cmd->args.collapse_polylines;
    470   tree_args.nthreads_hint = cmd->args.nthreads_hint;
    471 
    472   if(cmd->args.output) {
    473     write_args.filename = cmd->args.output;
    474   } else {
    475     write_args.filename = "stdout";
    476     write_args.file = stdout;
    477   }
    478 
    479   if((res = setup_tree_mixture(cmd, tree_args.molecules)) != RES_OK) goto error;
    480   if((res = sln_tree_create(cmd->sln, &tree_args, &tree)) != RES_OK) goto error;
    481   if((res = sln_tree_write(tree, &write_args)) != RES_OK) goto error;
    482 
    483 exit:
    484   if(tree) SLN(tree_ref_put(tree));
    485   return res;
    486 error:
    487   goto exit;
    488 }
    489 
    490 /*******************************************************************************
    491  * The program
    492  ******************************************************************************/
    493 int
    494 main(int argc, char** argv)
    495 {
    496   struct args args = ARGS_DEFAULT;
    497   struct cmd cmd = CMD_NULL;
    498   int err = 0;
    499   res_T res = RES_OK;
    500 
    501   if((res = args_init(&args, argc, argv)) != RES_OK) goto error;
    502   if(args.quit) goto exit;
    503 
    504   if((res = cmd_init(&cmd, &args)) != RES_OK) goto error;
    505   if((res = cmd_run(&cmd)) != RES_OK) goto error;
    506 
    507 exit:
    508   cmd_release(&cmd);
    509   CHK(mem_allocated_size() == 0);
    510   return err;
    511 error:
    512   err = 1;
    513   goto exit;
    514 }