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


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