star-line

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

commit 36625f9007d251b49309db36b22b09c5ca61cc42
parent fa18c937930ec2e14dce2d31e81cfee4be8310cf
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon, 20 Apr 2026 18:20:24 +0200

Add an accessor for the tree's spectral lines

Remove the accessor for a node's spectral lines. This function was a
double mistake. On the one hand, the spectral lines returned were those
from the spectroscopic database, not the spectral lines corresponding to
the thermodynamic properties for which the tree had been constructed. On
the other hand, there was no reason to query the line from the node. A
node encompasses a set of lines that, in any case, must be accessible
from the tree itself.

Consequently, the new function sln_tree_get_line returns an sln_line
based on its index in the spectroscopic database. This public structure
is actually the internal structure that this commit moved into the API.
The line properties are those corresponding to the thermodynamic
properties for which the tree was constructed.

A node's descriptor now contains the range of line indices it
encompasses. The number of lines has been removed from the descriptor
because it can be calculated from this range. The function returning the
number of lines in the node has also been removed because, again, this
number can be calculated from the node descriptor.

The tests and utilities have been updated to reflect these API changes.
Note that a new test has been added to verify that the value of a node
at a given wave number matches the sum of the contributions from its
lines. The new public function 'sln_line_eval' allows the caller to
evaluate the contribution of a line for a given wave number.

Diffstat:
Msrc/sln.h | 60+++++++++++++++++++++++++++++++++++++++---------------------
Msrc/sln_get.c | 4+++-
Msrc/sln_line.c | 71+++++++++++++++++++++++++++++++++++++----------------------------------
Msrc/sln_line.h | 18+-----------------
Msrc/sln_tree.c | 57++++++++++++++++++++++++++++++++-------------------------
Msrc/test_sln_tree.c | 115++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------
6 files changed, 191 insertions(+), 134 deletions(-)

diff --git a/src/sln.h b/src/sln.h @@ -196,11 +196,12 @@ struct sln_tree_desc { static const struct sln_tree_desc SLN_TREE_DESC_NULL = SLN_TREE_DESC_NULL__; struct sln_node_desc { - size_t nlines; + /* Range of lines belonging to the node. The endpoints are included */ + size_t ilines[2]; size_t nvertices; unsigned nchildren; }; -#define SLN_NODE_DESC_NULL__ {0,0,0} +#define SLN_NODE_DESC_NULL__ {{0,0},0,0} static const struct sln_node_desc SLN_NODE_DESC_NULL = SLN_NODE_DESC_NULL__; struct sln_vertex { /* 8 Bytes */ @@ -228,6 +229,16 @@ struct sln_mixture_load_args { static const struct sln_mixture_load_args SLN_MIXTURE_LOAD_ARGS_NULL = SLN_MIXTURE_LOAD_ARGS_NULL__; +struct sln_line { + double wavenumber; /* Line center wrt pressure in cm^-1 */ + double profile_factor; /* m^-1.cm^-1 (1e2*density*intensity) */ + double gamma_d; /* Doppler half width */ + double gamma_l; /* Lorentz half width */ + enum shtr_molecule_id molecule_id; +}; +#define SLN_LINE_NULL__ {0,0,0,0,SHTR_MOLECULE_ID_NULL} +static const struct sln_line SLN_LINE_NULL = SLN_LINE_NULL__; + /* Forward declarations of opaque data structures */ struct sln_device; struct sln_mixture; @@ -318,6 +329,20 @@ SLN_API const struct sln_node* /* NULL <=> No node */ sln_tree_get_root (const struct sln_tree* tree); +SLN_API res_T +sln_tree_get_line + (const struct sln_tree* tree, + const size_t iline, + struct sln_line* line); + +SLN_API res_T +sln_tree_write + (const struct sln_tree* tree, + const struct sln_tree_write_args* args); + +/******************************************************************************* + * Node API + ******************************************************************************/ SLN_API int sln_node_is_leaf (const struct sln_node* node); @@ -334,16 +359,17 @@ sln_node_get_child const struct sln_node* node, const unsigned ichild); /* 0 or #children */ -SLN_API size_t -sln_node_get_line_count - (const struct sln_node* node); +SLN_API double +sln_node_eval + (const struct sln_tree* tree, + const struct sln_node* node, + const double wavenumber); /* In cm^-1 */ SLN_API res_T -sln_node_get_line +sln_node_get_desc (const struct sln_tree* tree, const struct sln_node* node, - const size_t iline, - struct shtr_line* line); + struct sln_node_desc* desc); SLN_API res_T sln_node_get_mesh @@ -351,28 +377,20 @@ sln_node_get_mesh const struct sln_node* node, struct sln_mesh* mesh); +/******************************************************************************* + * Miscellaneous + ******************************************************************************/ SLN_API double -sln_node_eval +sln_line_eval (const struct sln_tree* tree, - const struct sln_node* node, + const struct sln_line* line, const double wavenumber); /* In cm^-1 */ -SLN_API res_T -sln_node_get_desc - (const struct sln_tree* tree, - const struct sln_node* node, - struct sln_node_desc* desc); - SLN_API double sln_mesh_eval (const struct sln_mesh* mesh, const double wavenumber); /* In cm^-1 */ -SLN_API res_T -sln_tree_write - (const struct sln_tree* tree, - const struct sln_tree_write_args* args); - /******************************************************************************* * Helper functions ******************************************************************************/ diff --git a/src/sln_get.c b/src/sln_get.c @@ -362,6 +362,7 @@ print_node_descriptor(const struct cmd* cmd) { const struct sln_node* node = NULL; struct sln_node_desc desc = SLN_NODE_DESC_NULL; + size_t nlines = 0; unsigned depth = 0; res_T res = RES_OK; ASSERT(cmd); @@ -371,8 +372,9 @@ print_node_descriptor(const struct cmd* cmd) res = sln_node_get_desc(cmd->tree, node, &desc); if(res != RES_OK) goto error; + nlines = desc.ilines[1] - desc.ilines[0] + 1/*inclusive bounds*/; printf("level: %u\n", depth); - printf("#lines: %lu\n", (unsigned long)desc.nlines); + printf("#lines: %lu\n", (unsigned long)nlines); printf("#vertices: %lu\n", (unsigned long)desc.nvertices); printf("#children: %u\n", desc.nchildren); diff --git a/src/sln_line.c b/src/sln_line.c @@ -184,7 +184,7 @@ error: static res_T regular_mesh_fragmented (const struct sln_tree* tree, - const struct line* line, + const struct sln_line* line, const size_t nvertices, struct darray_double* wavenumbers) /* List of issued vertices */ { @@ -248,7 +248,7 @@ error: static res_T eval_mesh (const struct sln_tree* tree, - const struct line* line, + const struct sln_line* line, const struct darray_double* wavenumbers, struct darray_double* values) { @@ -267,7 +267,7 @@ eval_mesh nu = darray_double_cdata_get(wavenumbers); ha = darray_double_data_get(values); FOR_EACH(ivertex, 0, nvertices) { - ha[ivertex] = line_value(tree, line, nu[ivertex]); + ha[ivertex] = sln_line_eval(tree, line, nu[ivertex]); } exit: @@ -346,7 +346,7 @@ next_vertex_value static res_T save_line_mesh (struct sln_tree* tree, - const struct line* line, + const struct sln_line* line, const struct darray_double* wavenumbers, const struct darray_double* values, struct darray_vertex* vertices, /* buffer in which vertices are added */ @@ -427,7 +427,7 @@ res_T line_setup (const struct sln_tree* tree, const size_t iline, - struct line* line) + struct sln_line* line) { struct shtr_molecule molecule = SHTR_MOLECULE_NULL; struct shtr_line shtr_line = SHTR_LINE_NULL; @@ -466,34 +466,6 @@ error: goto exit; } -double -line_value - (const struct sln_tree* tree, - const struct line* line, - const double wavenumber) -{ - const struct sln_molecule* mol_params = NULL; - double profile = 0; - ASSERT(tree && line); - - /* Retrieve the molecular parameters of the line to be mesh */ - mol_params = tree->args.molecules + line->molecule_id; - - if(wavenumber < line->wavenumber - mol_params->cutoff - || wavenumber > line->wavenumber + mol_params->cutoff) { - return 0; - } - - switch(tree->args.line_profile) { - case SLN_LINE_PROFILE_VOIGT: - profile = sln_compute_voigt_profile - (wavenumber, line->wavenumber, line->gamma_d, line->gamma_l); - break; - default: FATAL("Unreachable code.\n"); break; - } - return line->profile_factor * profile; -} - res_T line_mesh (struct sln_tree* tree, @@ -503,7 +475,7 @@ line_mesh size_t vertices_range[2]) /* out. Bounds are inclusive */ { /* The line */ - struct line line = LINE_NULL; + struct sln_line line = SLN_LINE_NULL; /* Temporary mesh */ struct darray_double values; /* List of evaluated values */ @@ -558,3 +530,34 @@ exit: error: goto exit; } + +/******************************************************************************* + * Exported functions + ******************************************************************************/ +double +sln_line_eval + (const struct sln_tree* tree, + const struct sln_line* line, + const double wavenumber) +{ + const struct sln_molecule* mol_params = NULL; + double profile = 0; + ASSERT(tree && line); + + /* Retrieve the molecular parameters of the line to be mesh */ + mol_params = tree->args.molecules + line->molecule_id; + + if(wavenumber < line->wavenumber - mol_params->cutoff + || wavenumber > line->wavenumber + mol_params->cutoff) { + return 0; + } + + switch(tree->args.line_profile) { + case SLN_LINE_PROFILE_VOIGT: + profile = sln_compute_voigt_profile + (wavenumber, line->wavenumber, line->gamma_d, line->gamma_l); + break; + default: FATAL("Unreachable code.\n"); break; + } + return line->profile_factor * profile; +} diff --git a/src/sln_line.h b/src/sln_line.h @@ -25,16 +25,6 @@ #include <rsys/rsys.h> #include <math.h> -struct line { - double wavenumber; /* Line center wrt pressure in cm^-1 */ - double profile_factor; /* m^-1.cm^-1 (1e2*density*intensity) */ - double gamma_d; /* Doppler half width */ - double gamma_l; /* Lorentz half width */ - enum shtr_molecule_id molecule_id; -}; -#define LINE_NULL__ {0} -static const struct line LINE_NULL = LINE_NULL__; - /* Forward declaration */ struct sln_mixture; struct sln_tree; @@ -53,13 +43,7 @@ extern LOCAL_SYM res_T line_setup (const struct sln_tree* tree, const size_t iline, - struct line* line); - -extern LOCAL_SYM double -line_value - (const struct sln_tree* tree, - const struct line* line, - const double wavenumber); + struct sln_line* line); extern LOCAL_SYM res_T line_mesh diff --git a/src/sln_tree.c b/src/sln_tree.c @@ -626,6 +626,33 @@ sln_tree_get_root(const struct sln_tree* tree) } } +res_T +sln_tree_get_line + (const struct sln_tree* tree, + const size_t iline, + struct sln_line* line) +{ + size_t nlines = 0; + res_T res = RES_OK; + + if(!tree || !line) { res = RES_BAD_ARG; goto error; } + + SHTR(line_list_get_size(tree->args.lines, &nlines)); + if(iline >= nlines) { res = RES_BAD_ARG; goto error; } + + res = line_setup(tree, iline, line); + if(res != RES_OK) { + ERROR(tree->sln, "%s: could not setup the line %lu-- %s\n", + FUNC_NAME, iline, res_to_cstr(res)); + goto error; + } + +exit: + return res; +error: + goto exit; +} + int sln_node_is_leaf(const struct sln_node* node) { @@ -659,27 +686,6 @@ sln_node_get_child return node + node->offset + ichild; } -size_t -sln_node_get_line_count(const struct sln_node* node) -{ - ASSERT(node); - return node->range[1] - node->range[0] + 1/*Both boundaries are inclusives*/; -} - -res_T -sln_node_get_line - (const struct sln_tree* tree, - const struct sln_node* node, - const size_t iline, - struct shtr_line* line) -{ - if(!tree || !node || iline > sln_node_get_line_count(node)) - return RES_BAD_ARG; - - return shtr_line_list_at - (tree->args.lines, node->range[0] + iline, line); -} - res_T sln_node_get_mesh (const struct sln_tree* tree, @@ -703,16 +709,17 @@ sln_node_eval ASSERT(tree && node); FOR_EACH(iline, node->range[0], node->range[1]+1) { - struct line line = LINE_NULL; + struct sln_line line = SLN_LINE_NULL; res_T res = RES_OK; res = line_setup(tree, iline, &line); if(res != RES_OK) { WARN(tree->sln, "%s: could not setup the line %lu-- %s\n", FUNC_NAME, iline, res_to_cstr(res)); + continue; } - ka += line_value(tree, &line, nu); + ka += sln_line_eval(tree, &line, nu); } return ka; } @@ -724,8 +731,8 @@ sln_node_get_desc struct sln_node_desc* desc) { if(!tree || !node || !desc) return RES_BAD_ARG; - desc->nlines = node->range[1] - node->range[0]; - desc->nlines += 1/*boundaries are inclusives*/; + desc->ilines[0] = node->range[0]; + desc->ilines[1] = node->range[1]; desc->nvertices = node->nvertices; desc->nchildren = sln_node_get_child_count(tree, node); return RES_OK; diff --git a/src/test_sln_tree.c b/src/test_sln_tree.c @@ -101,11 +101,17 @@ check_tree_lines node = sln_tree_get_root(tree); while(node) { + struct sln_node_desc node_desc = SLN_NODE_DESC_NULL; + size_t node_nlines = 0; + + CHK(sln_node_get_desc(tree, node, &node_desc) == RES_OK); + node_nlines = node_desc.ilines[1] - node_desc.ilines[0] + 1; + if(!sln_node_is_leaf(node)) { unsigned i = 0; unsigned n = sln_node_get_child_count(tree, node); - CHK(sln_node_get_line_count(node) > desc.leaf_nlines); + CHK(node_nlines > desc.leaf_nlines); FOR_EACH(i, 1, n) { stack[istack++] = sln_node_get_child(tree, node, i); @@ -113,20 +119,13 @@ check_tree_lines node = sln_node_get_child(tree, node, 0); /* Handle the child 0 */ } else { - size_t iline, nlines; - - nlines = sln_node_get_line_count(node); - CHK(nlines <= desc.leaf_nlines); - FOR_EACH(iline, 0, nlines) { - struct shtr_line line = SHTR_LINE_NULL; - size_t found_iline = 0; - CHK(sln_node_get_line(tree, node, iline, &line) == RES_OK); - found_iline = find_line(line_list, &line); - CHK(found_iline != SIZE_MAX); /* Line is found */ - - if(!found_lines[found_iline]) { + size_t iline; + + CHK(node_nlines <= desc.leaf_nlines); + FOR_EACH(iline, node_desc.ilines[0], node_desc.ilines[1]+1) { + if(!found_lines[iline]) { found_nlines += 1; - found_lines[found_iline] = 1; + found_lines[iline] = 1; } } @@ -163,19 +162,31 @@ check_node_equality { struct sln_mesh mesh1 = SLN_MESH_NULL; struct sln_mesh mesh2 = SLN_MESH_NULL; - size_t iline, nlines; + struct sln_node_desc desc1 = SLN_NODE_DESC_NULL; + struct sln_node_desc desc2 = SLN_NODE_DESC_NULL; + size_t iline = 0; CHK(node1 && node2); CHK(sln_node_is_leaf(node1) == sln_node_is_leaf(node2)); - CHK(sln_node_get_line_count(node1) == sln_node_get_line_count(node2)); - nlines = sln_node_get_line_count(node1); - - FOR_EACH(iline, 0, nlines) { - struct shtr_line line1 = SHTR_LINE_NULL; - struct shtr_line line2 = SHTR_LINE_NULL; - CHK(sln_node_get_line(tree1, node1, iline, &line1) == RES_OK); - CHK(sln_node_get_line(tree2, node2, iline, &line2) == RES_OK); - CHK(shtr_line_eq(&line1, &line2)); + CHK(sln_node_get_desc(tree1, node1, &desc1) == RES_OK); + CHK(sln_node_get_desc(tree2, node2, &desc2) == RES_OK); + + CHK(desc1.ilines[0] == desc2.ilines[0]); + CHK(desc1.ilines[1] == desc2.ilines[1]); + CHK(desc1.nvertices == desc2.nvertices); + CHK(desc1.nchildren == desc2.nchildren); + + FOR_EACH(iline, desc1.ilines[0], desc1.ilines[1]+1) { + struct sln_line line1 = SLN_LINE_NULL; + struct sln_line line2 = SLN_LINE_NULL; + CHK(sln_tree_get_line(tree1, iline, &line1) == RES_OK); + CHK(sln_tree_get_line(tree2, iline, &line2) == RES_OK); + + CHK(line1.wavenumber == line2.wavenumber); + CHK(line1.profile_factor == line2.profile_factor); + CHK(line1.gamma_d == line2.gamma_d); + CHK(line1.gamma_l == line2.gamma_l); + CHK(line1.molecule_id == line2.molecule_id); } CHK(sln_node_get_mesh(tree1, node1, &mesh1) == RES_OK); @@ -245,6 +256,33 @@ check_tree_equality } static void +check_node_value(const struct sln_tree* tree, const struct sln_node* node) +{ + struct sln_node_desc desc = SLN_NODE_DESC_NULL; + struct sln_line line = SLN_LINE_NULL; + size_t iline = 0; + double ka_ref = 0; + double ka_node = 0; + double nu = 0; + + CHK(sln_node_get_desc(tree, node, &desc) == RES_OK); + + CHK(sln_tree_get_line(tree, desc.ilines[0], &line) == RES_OK); + nu = line.wavenumber + 10; + CHK(sln_tree_get_line(tree, desc.ilines[1], &line) == RES_OK); + nu += line.wavenumber - 10; + nu *= 0.5; + + ka_node = sln_node_eval(tree, node, nu); + FOR_EACH(iline, desc.ilines[0], desc.ilines[1]+1/*inclusive*/) { + CHK(sln_tree_get_line(tree, iline, &line) == RES_OK); + ka_ref += sln_line_eval(tree, &line, nu); + } + + CHK(eq_eps(ka_node, ka_ref, ka_ref*1e-6)); +} + +static void test_tree (struct sln_device* sln, const struct sln_tree_create_args* tree_args_in, @@ -256,7 +294,7 @@ test_tree struct sln_mesh mesh = SLN_MESH_NULL; struct sln_tree* tree = NULL; const struct sln_node* node = NULL; - struct shtr_line line = SHTR_LINE_NULL; + struct sln_line line = SLN_LINE_NULL; size_t nlines = 0; size_t nleaf = 0; unsigned depth = 0; @@ -291,6 +329,12 @@ test_tree CHK(desc.temperature == tree_args.temperature); CHK(desc.arity == tree_args.arity); + CHK(sln_tree_get_line(NULL, 0, &line) == RES_BAD_ARG); + CHK(sln_tree_get_line(tree, nlines, &line) == RES_BAD_ARG); + CHK(sln_tree_get_line(tree, 0, NULL) == RES_BAD_ARG); + CHK(sln_tree_get_line(tree, 0, &line) == RES_OK); + CHK(sln_tree_get_line(tree, nlines-1, &line) == RES_OK); + CHK(node = sln_tree_get_root(tree)); CHK(node != NULL); @@ -298,7 +342,8 @@ test_tree CHK(sln_node_get_desc(tree, NULL, &node_desc) == RES_BAD_ARG); CHK(sln_node_get_desc(tree, node, NULL) == RES_BAD_ARG); CHK(sln_node_get_desc(tree, node, &node_desc) == RES_OK); - CHK(node_desc.nlines = desc.nlines); + CHK(node_desc.ilines[0] == 0); + CHK(node_desc.ilines[1] == desc.nlines - 1); CHK(node_desc.nvertices >= 1 && node_desc.nvertices < desc.nvertices); CHK(node_desc.nchildren <= desc.arity); @@ -311,8 +356,12 @@ test_tree node = sln_node_get_child(tree, node, 0); CHK(sln_node_get_desc(tree, node, &node_desc_next) == RES_OK); - CHK(node_desc_next.nlines >= 1); - CHK(node_desc_next.nlines < node_desc.nlines); + + #define NLINES(Desc) ((Desc).ilines[1] - (Desc).ilines[0] + 1) + CHK(NLINES(node_desc_next) >= 1); + CHK(NLINES(node_desc_next) < NLINES(node_desc)); + #undef NLINES + CHK(node_desc_next.nvertices >= 1); if(sln_node_is_leaf(node)) { CHK(node_desc_next.nchildren == 0); @@ -325,19 +374,13 @@ test_tree } CHK(sln_node_get_child_count(tree, node) == 0); - CHK(sln_node_get_line_count(node) <= desc.leaf_nlines); - CHK(sln_node_get_line(NULL, node, 0, &line) == RES_BAD_ARG); - CHK(sln_node_get_line(tree, NULL, 0, &line) == RES_BAD_ARG); - CHK(sln_node_get_line(tree, node, desc.leaf_nlines+1, &line) - == RES_BAD_ARG); - CHK(sln_node_get_line(tree, node, 0, NULL) == RES_BAD_ARG); - CHK(sln_node_get_line(tree, node, 0, &line) == RES_OK); - CHK(sln_node_get_line(tree, sln_tree_get_root(tree), 0, &line) == RES_OK); + CHK(node_desc.ilines[1] - node_desc.ilines[0] + 1 <= desc.leaf_nlines); CHK(sln_node_get_mesh(NULL, node, &mesh) == RES_BAD_ARG); CHK(sln_node_get_mesh(tree, NULL, &mesh) == RES_BAD_ARG); CHK(sln_node_get_mesh(tree, node, NULL) == RES_BAD_ARG); CHK(sln_node_get_mesh(tree, node, &mesh) == RES_OK); + check_node_value(tree, node); CHK(node = sln_tree_get_root(tree)); check_tree_lines(tree, line_list);