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 (13307B)


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