star-line

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

commit a501a85b20aeb2bc70475664ec4f66597f7c7ab9
parent 16bc9476d789578cacdb255cd9fa736a48198c06
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 23 Sep 2026 11:14:23 +0200

Test for different thermo conditions

This new test checks that although built under specific
thermodynamic conditions, the caller can calculate the value of a tree
node and a line which it partitions under different conditions.

The same tree can therefore be used to sample the lines by importance at
other thermodynamic conditions as long as the criteria for such sampling
are met, such as the node polylines which must be an upper bound of the
lines it encompasses.

Diffstat:
MMakefile | 2++
Asrc/test_sln_thermo_props.c | 377+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 files changed, 379 insertions(+), 0 deletions(-)

diff --git a/Makefile b/Makefile @@ -212,6 +212,7 @@ TEST_SRC =\ src/test_sln_device.c\ src/test_sln_mesh.c\ src/test_sln_mixture.c\ + src/test_sln_thermo_props.c\ src/test_sln_tree.c\ src/test_sln_tree_sample.c TEST_OBJ = $(TEST_SRC:.c=.o) @@ -246,6 +247,7 @@ $(TEST_OBJ): config.mk sln-local.pc test_sln_device \ test_sln_mesh \ test_sln_mixture \ +test_sln_thermo_props\ test_sln_tree_sample \ : config.mk sln-local.pc $(LIBNAME) $(CC) $(CFLAGS_TEST) -o $@ src/$@.o $(LDFLAGS_TEST) diff --git a/src/test_sln_thermo_props.c b/src/test_sln_thermo_props.c @@ -0,0 +1,377 @@ +/* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com) + * Copyright (C) 2026 Université de Lorraine + * Copyright (C) 2022 Centre National de la Recherche Scientifique + * Copyright (C) 2022 Université Paul Sabatier + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#include <test_sln_lines.h> + +#include "sln.h" + +#include <rsys/mem_allocator.h> +#include <rsys/rsys.h> + +struct thermo_props { + double xH2O; + double xCO2; + double xO3; + double pressure; /*[atm]*/ + double temperature; /*[K]*/ +}; + +static const struct thermo_props thermo_props1 = {0.15, 0.10, 0.05, 10, 600}; +static const struct thermo_props thermo_props2 = {0.10, 0.15, 0.25, 20, 450}; + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static struct shtr_isotope_metadata* +setup_isotopes(struct shtr* shtr) +{ + struct shtr_isotope_metadata* metadata = NULL; + FILE* fp = NULL; + + CHK(fp = tmpfile()); + fprintf(fp, "Molecule # Iso Abundance Q(296K) gj Molar Mass(g)\n"); + write_shtr_molecule(fp, &g_H2O); + write_shtr_molecule(fp, &g_CO2); + write_shtr_molecule(fp, &g_O3); + rewind(fp); + + CHK(shtr_isotope_metadata_load_stream(shtr, fp, NULL, &metadata) == RES_OK); + + CHK(fclose(fp) == 0); + + return metadata; +} + +static struct shtr_line_list* +setup_lines(struct shtr* shtr) +{ + struct shtr_line_list_load_args args = SHTR_LINE_LIST_LOAD_ARGS_NULL; + struct shtr_line_list* lines = NULL; + FILE* fp = NULL; + + CHK(fp = tmpfile()); + write_shtr_lines(fp, g_lines, g_nlines); + rewind(fp); + + args.filename = "stream"; + args.file = fp; + CHK(shtr_line_list_load(shtr, &args, &lines) == RES_OK); + + CHK(fclose(fp) == 0); + + return lines; +} + +static struct sln_tree* +create_tree + (struct sln_device* sln, + struct shtr_isotope_metadata* mdata, + struct shtr_line_list* lines, + const struct thermo_props* props) +{ + struct sln_tree_create_args args = SLN_TREE_CREATE_ARGS_DEFAULT; + struct sln_tree* tree = NULL; + + args.metadata = mdata; + args.lines = lines; + + args.molecules[SHTR_H2O].concentration = props->xH2O; + args.molecules[SHTR_H2O].cutoff = 25; /* [cm^-1] */ + args.molecules[SHTR_CO2].concentration = props->xCO2; + args.molecules[SHTR_CO2].cutoff = 50; /* [cm^-1] */ + args.molecules[SHTR_O3].concentration = props->xO3; + args.molecules[SHTR_O3].cutoff = 25; /* [cm^-1] */ + + args.pressure = props->pressure; /*[atm]*/ + args.temperature = props->temperature; /*[K]*/ + + CHK(sln_tree_create(sln, &args, &tree) == RES_OK); + return tree; +} + +static INLINE double /* in [0,1[ */ +rand_canonic(void) +{ + return (double)rand() / (double)((long)RAND_MAX+1); +} + +static double /* [cm^-1] */ +line_sample_nu(const struct sln_tree* tree, const size_t iline) +{ + struct sln_line line = SLN_LINE_NULL; + double nu_range[2] = {0,0}; /* [cm^-1] */ + double nu = 0; /* [cm^-1] */ + + CHK(sln_tree_get_line(tree, iline, NULL, &line) == RES_OK); + nu_range[0] = line.wavenumber - 20; + nu_range[1] = line.wavenumber + 20; + nu = nu_range[0] + rand_canonic() * (nu_range[1] - nu_range[0]); + return nu; +} + +static double /* [cm^-1] */ +node_sample_nu(const struct sln_tree* tree, const struct sln_node* node) +{ + struct sln_line line = SLN_LINE_NULL; + struct sln_node_desc desc = SLN_NODE_DESC_NULL; + double nu_range[2] = {0,0}; /* [cm^-1] */ + double nu = 0; /* [cm^-1] */ + + CHK(sln_node_get_desc(tree, node, &desc) == RES_OK); + + CHK(sln_tree_get_line(tree, desc.ilines[0], NULL, &line) == RES_OK); + nu_range[0] = line.wavenumber - 10; + CHK(sln_tree_get_line(tree, desc.ilines[1], NULL, &line) == RES_OK); + nu_range[1] = line.wavenumber + 10; + + nu = nu_range[0] + rand_canonic() * (nu_range[1] - nu_range[0]); + return nu; +} + +static size_t +node_sample_line(const struct sln_tree* tree, const struct sln_node* node) +{ + struct sln_node_desc desc = SLN_NODE_DESC_NULL; + const double r = rand_canonic(); + size_t iline = 0; + + CHK(sln_node_get_desc(tree, node, &desc) == RES_OK); + iline = desc.ilines[0]; + iline += (size_t)(r * (double)(desc.ilines[1] - desc.ilines[0] + 1)); + return iline; +} + +/* Check that, even belonging to 2 trees built from different thermodynamic + * properties, a line has the _exact_ same value when queried with the same + * thermodynamic properties. */ +static void +cmp_lines_values + (const struct sln_tree* tree1, + const struct sln_tree* tree2, + const size_t iline, + const double nu/*[cm^-1]*/) +{ + struct sln_line line1 = SLN_LINE_NULL; + struct sln_line line2 = SLN_LINE_NULL; + struct sln_thermo_props props = SLN_THERMO_PROPS_NULL; + double ka1 = 0; + double ka2 = 0; + + CHK(sln_tree_get_line(tree1, iline, NULL, &line1) == RES_OK); + CHK(sln_tree_get_line(tree2, iline, NULL, &line2) == RES_OK); + ka1 = sln_line_eval(tree1, &line1, nu); + ka2 = sln_line_eval(tree2, &line2, nu); + CHK(ka1 != ka2); + + props.concentrations[SHTR_H2O] = thermo_props1.xH2O; + props.concentrations[SHTR_CO2] = thermo_props1.xCO2; + props.concentrations[SHTR_O3] = thermo_props1.xO3; + props.pressure = thermo_props1.pressure; + props.temperature = thermo_props1.temperature; + CHK(sln_tree_get_line(tree1, iline, NULL, &line1) == RES_OK); + CHK(sln_tree_get_line(tree2, iline, &props, &line2) == RES_OK); + ka1 = sln_line_eval(tree1, &line1, nu); + ka2 = sln_line_eval(tree2, &line2, nu); + CHK(ka1 == ka2); + + props.concentrations[SHTR_H2O] = thermo_props2.xH2O; + props.concentrations[SHTR_CO2] = thermo_props2.xCO2; + props.concentrations[SHTR_O3] = thermo_props2.xO3; + props.pressure = thermo_props2.pressure; + props.temperature = thermo_props2.temperature; + CHK(sln_tree_get_line(tree1, iline, &props, &line1) == RES_OK); + CHK(sln_tree_get_line(tree2, iline, NULL, &line2) == RES_OK); + ka1 = sln_line_eval(tree1, &line1, nu); + ka2 = sln_line_eval(tree2, &line2, nu); + CHK(ka1 == ka2); + + props.concentrations[SHTR_H2O] = thermo_props1.xH2O; + props.concentrations[SHTR_CO2] = thermo_props2.xCO2; + props.concentrations[SHTR_O3] = thermo_props1.xO3; + props.pressure = thermo_props2.pressure; + props.temperature = thermo_props1.temperature; + CHK(sln_tree_get_line(tree1, iline, &props, &line1) == RES_OK); + CHK(sln_tree_get_line(tree2, iline, &props, &line2) == RES_OK); + ka1 = sln_line_eval(tree1, &line1, nu); + ka2 = sln_line_eval(tree2, &line2, nu); + CHK(ka1 == ka2); +} + +static void +cmp_lines + (const struct sln_tree* tree1, + const struct sln_tree* tree2, + const size_t iline) +{ + const size_t N = 50; + size_t i = 0; + + FOR_EACH(i, 0, N) { + const double nu = line_sample_nu(tree1, iline); + cmp_lines_values(tree1, tree2, iline, nu); + } +} + +static void +cmp_nodes_lines + (const struct sln_tree* tree1, + const struct sln_node* node1, + const struct sln_tree* tree2, + const struct sln_node* node2) +{ + const size_t N = 50; + size_t i = 0; + (void)node2; + + FOR_EACH(i, 0, N) { + const size_t iline = node_sample_line(tree1, node1); + cmp_lines(tree1, tree2, iline); + } +} + +/* Check that, even belonging to 2 trees built from different thermodynamic + * properties, a node has the _exact_ same value when evaluated with the same + * thermodynamic properties. */ +static void +cmp_nodes_values + (const struct sln_tree* tree1, + const struct sln_node* node1, + const struct sln_tree* tree2, + const struct sln_node* node2, + const double nu/*[cm^-1]*/) +{ + struct sln_thermo_props props = SLN_THERMO_PROPS_NULL; + double ka1 = 0; + double ka2 = 0; + + ka1 = sln_node_eval(tree1, node1, NULL, nu); + ka2 = sln_node_eval(tree2, node2, NULL, nu); + CHK(ka1 != ka2); + + props.concentrations[SHTR_H2O] = thermo_props1.xH2O; + props.concentrations[SHTR_CO2] = thermo_props1.xCO2; + props.concentrations[SHTR_O3] = thermo_props1.xO3; + props.pressure = thermo_props1.pressure; + props.temperature = thermo_props1.temperature; + ka1 = sln_node_eval(tree1, node1, NULL, nu); + ka2 = sln_node_eval(tree2, node2, &props, nu); + CHK(ka1 == ka2); + + props.concentrations[SHTR_H2O] = thermo_props2.xH2O; + props.concentrations[SHTR_CO2] = thermo_props2.xCO2; + props.concentrations[SHTR_O3] = thermo_props2.xO3; + props.pressure = thermo_props2.pressure; + props.temperature = thermo_props2.temperature; + ka1 = sln_node_eval(tree1, node1, &props, nu); + ka2 = sln_node_eval(tree2, node2, NULL, nu); + CHK(ka1 == ka2); + + props.concentrations[SHTR_H2O] = thermo_props1.xH2O; + props.concentrations[SHTR_CO2] = thermo_props2.xCO2; + props.concentrations[SHTR_O3] = thermo_props1.xO3; + props.pressure = thermo_props2.pressure; + props.temperature = thermo_props1.temperature; + ka1 = sln_node_eval(tree1, node1, &props, nu); + ka2 = sln_node_eval(tree2, node2, &props, nu); + CHK(ka1 == ka2); +} + +static void +cmp_nodes + (const struct sln_tree* tree1, + const struct sln_node* node1, + const struct sln_tree* tree2, + const struct sln_node* node2) +{ + const size_t N = 50; + size_t i = 0; + + FOR_EACH(i, 0, N) { + const double nu = node_sample_nu(tree1, node1); + cmp_nodes_values(tree1, node1, tree2, node2, nu); + } + + cmp_nodes_lines(tree1, node1, tree2, node2); +} + +static void +cmp_trees(const struct sln_tree* tree1, const struct sln_tree* tree2) +{ + const struct sln_node* node1 = NULL; + const struct sln_node* node2 = NULL; + struct sln_tree_desc desc = SLN_TREE_DESC_NULL; + + node1 = sln_tree_get_root(tree1); + node2 = sln_tree_get_root(tree2); + + CHK(sln_tree_get_desc(tree1, &desc) == RES_OK); + CHK(desc.arity == 2); /* Assume that the arity of the tree is 2 */ + + for(;;) { + unsigned ichild = 0; + + cmp_nodes(tree1, node1, tree2, node2); + if(sln_node_is_leaf(node1)) break; + + /* Randomly choose one node child */ + ichild = rand_canonic() < 0.5 ? 0 : 1; + node1 = sln_node_get_child(tree1, node1, ichild); + node2 = sln_node_get_child(tree2, node2, ichild); + } +} + +/******************************************************************************* + * The test + ******************************************************************************/ +int +main(void) +{ + struct shtr_create_args shtr_args = SHTR_CREATE_ARGS_DEFAULT; + struct shtr* shtr = NULL; + struct shtr_isotope_metadata* mdata = NULL; + struct shtr_line_list* lines = NULL; + + struct sln_device_create_args sln_args = SLN_DEVICE_CREATE_ARGS_DEFAULT; + struct sln_device* sln = NULL; + struct sln_tree* tree1 = NULL; + struct sln_tree* tree2 = NULL; + + shtr_args.verbose = 1; + CHK(shtr_create(&shtr_args, &shtr) == RES_OK); + sln_args.verbose = 1; + CHK(sln_device_create(&sln_args, &sln) == RES_OK); + + mdata = setup_isotopes(shtr); + lines = setup_lines(shtr); + + tree1 = create_tree(sln, mdata, lines, &thermo_props1); + tree2 = create_tree(sln, mdata, lines, &thermo_props2); + + cmp_trees(tree1, tree2); + + CHK(shtr_ref_put(shtr) == RES_OK); + CHK(shtr_line_list_ref_put(lines) == RES_OK); + CHK(shtr_isotope_metadata_ref_put(mdata) == RES_OK); + + CHK(sln_device_ref_put(sln) == RES_OK); + CHK(sln_tree_ref_put(tree1) == RES_OK); + CHK(sln_tree_ref_put(tree2) == RES_OK); + + CHK(mem_allocated_size() == 0); + return 0; +}