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


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