star-line

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

commit e667ad29cad1d9d9624a76a989671ce40d51a86a
parent a95c0a13fa08100f28f1f6ef15bc01fb83d2f33b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu, 29 Jan 2026 08:53:05 +0100

Rm the spectral range as an input argument for tree creation

Once again (see 592c934), the number of lines can be enormous, so it is
not wise to load them all before selecting a subset by limiting the
spectral range. The lines to be partitioned should only be defined from
the list of lines submitted as input. The spectral range of the tree
then corresponds to that of the lines, without the user having to define
it.

This commit also removes the constraint that wave numbers must be
positive. Now that the spectral range is no longer limited, lines can
extend to negative wave numbers. This is not a problem. First, because
it does not introduce any errors, since the user can always reduce the
interval they use, for example to only the positive values of the
spectrum. Furthermore, if some vertices are added unnecessarily, the
meshing of the lines is much simpler, with symmetric meshing in all
situations. Thus, a strategy aimed at reducing memory footprint could
use this symmetry, or simply store half of the line vertices. In other
words, the limited increase in memory usage introduced by this
commit allows for the implementation of strategies aimed at reducing
this same memory usage that are disproportionate to this small
additional cost.

Diffstat:
Msrc/sln.h | 3---
Msrc/sln_line.c | 129++++++++-----------------------------------------------------------------------
Msrc/sln_line.h | 2+-
Msrc/sln_tree.c | 6++----
Msrc/sln_tree_build.c | 2+-
Msrc/test_sln_tree.c | 2--
6 files changed, 17 insertions(+), 127 deletions(-)

diff --git a/src/sln.h b/src/sln.h @@ -103,8 +103,6 @@ struct sln_tree_create_args { * coarser the mesh */ float mesh_decimation_err; /* > 0 */ enum sln_mesh_type mesh_type; /* Type of mesh to generate */ - - double wavenumber_range[2]; /* [cm^-1, cm^-1] */ }; #define SLN_TREE_CREATE_ARGS_DEFAULT__ { \ NULL, /* metadata */ \ @@ -116,7 +114,6 @@ struct sln_tree_create_args { 16, /* #vertices hint */ \ 0.01f, /* Mesh decimation error */ \ SLN_MESH_UPPER, /* Mesh type */ \ - {-DBL_MAX, DBL_MAX} /* Spectral range */ \ } static const struct sln_tree_create_args SLN_TREE_CREATE_ARGS_DEFAULT = SLN_TREE_CREATE_ARGS_DEFAULT__; diff --git a/src/sln_line.c b/src/sln_line.c @@ -251,8 +251,6 @@ eval_mesh const double* nu = NULL; double* ha = NULL; size_t ivertex, nvertices; - double range[2]; - double cutoff; res_T res = RES_OK; ASSERT(tree && line && wavenumbers && values); @@ -262,34 +260,10 @@ eval_mesh res = darray_double_resize(values, nvertices); if(res != RES_OK) goto error; - /* Calculate the spectral range in which the vertices should be evaluated. - * Recall that a line being symmetrical in its center, we will only evaluate - * its upper half. Anyway, note that a line clipped by the spectral range is - * no longer symmetric. In the following, we thus adjust the spectral range - * to emit all the vertices needed to mesh the upper half _and_ the lower - * half of a possibly clipped line */ - cutoff = MMAX - (line->wavenumber - tree->args.wavenumber_range[0], - tree->args.wavenumber_range[1] - line->wavenumber); - range[0] = MMIN(line->wavenumber, tree->args.wavenumber_range[0]); - range[1] = line->wavenumber + cutoff; - nu = darray_double_cdata_get(wavenumbers); ha = darray_double_data_get(values); FOR_EACH(ivertex, 0, nvertices) { - const double nu_curr = nu[ivertex]; - const double nu_next = ivertex < nvertices - 1 ? nu[ivertex+1] : nu_curr; - const double nu_prev = ivertex > 0 ? nu[ivertex-1] : nu_curr; - - /* Check if the vertex needs to be cut. Note that we do not evaluate a vertex - * if it is out of spectral range and does not directly precede or succeed - * a vertex that is not clipped. This ensures that the spectral boundaries - * are always surrounded by vertices whose value has been evaluate. */ - if(nu_next < range[0] || range[1] < nu_prev) { - ha[ivertex] = -1; /* Clip */ - } else { - ha[ivertex] = line_value(tree, line, nu_curr); - } + ha[ivertex] = line_value(tree, line, nu[ivertex]); } exit: @@ -381,15 +355,13 @@ save_line_mesh const struct darray_double* values, size_t vertices_range[2]) /* Range into which the line vertices are saved */ { - const struct sln_molecule* mol_params = NULL; const double* wnums = NULL; const double* vals = NULL; - const double* spectral_range = NULL; - size_t nvertices; - size_t nwavenumbers; - size_t line_nvertices; - size_t ivertex; - size_t i; + size_t nvertices = 0; + size_t nwavenumbers = 0; + size_t line_nvertices = 0; + size_t ivertex = 0; + size_t i = 0; res_T res = RES_OK; ASSERT(tree && line && wavenumbers && values && vertices_range); @@ -397,7 +369,6 @@ save_line_mesh nvertices = darray_vertex_size_get(&tree->vertices); nwavenumbers = darray_double_size_get(wavenumbers); - spectral_range = tree->args.wavenumber_range; /* Compute the overall number of vertices of the line */ line_nvertices = nwavenumbers @@ -410,112 +381,38 @@ save_line_mesh wnums = darray_double_cdata_get(wavenumbers); vals = darray_double_cdata_get(values); - i = nvertices; - mol_params = tree->args.molecules + line->molecule_id; + i = nvertices; #define MIRROR(Nu) (2*line->wavenumber - (Nu)) - /* Emit the 1st vertex for the lower bound of the spectral range if the line - * is clipped by it */ - if(line->wavenumber - mol_params->cutoff <= spectral_range[0]) { - struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices)+i; - const double nu = spectral_range[0]; - double ha = 0; - - switch(tree->args.mesh_type) { - case SLN_MESH_UPPER: - ha = nu > line->wavenumber - ? next_vertex_value(nu, wavenumbers, values) - : next_vertex_value(MIRROR(nu), wavenumbers, values); - break; - - case SLN_MESH_USUAL: - ha = line_value(tree, line, nu); - break; - - default: FATAL("Unreachable code.\n"); break; - } - vtx->wavenumber = (float)nu; - vtx->ka = (float)ha; - ++i; - } - /* Copy the vertices of the line for its lower half */ FOR_EACH_REVERSE(ivertex, nwavenumbers-1, 0) { - struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices)+i; + struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices) + i++; const double nu = MIRROR(wnums[ivertex]); const double ha = vals[ivertex]; - /* The wavenumber is less than the lower bound of the spectral range. - * Discard the vertex */ - if(nu <= spectral_range[0]) continue; - - /* The wavenumber is greater than the upper bound of the spectral range. - * Stop the copy since the remaining vertices are greater than the current - * one and are thus necessary out of the spectral range */ - if(nu >= spectral_range[1]) break; - vtx->wavenumber = (float)nu; vtx->ka = (float)ha; - ++i; } /* Copy the vertices of the line for its upper half */ FOR_EACH(ivertex, 0, nwavenumbers) { - struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices)+i; + struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices) + i++; const double nu = wnums[ivertex]; const double ha = vals[ivertex]; - /* The wavenumber is less than the lower bound of the spectral range. - * Discard the vertex */ - if(nu <= spectral_range[0]) continue; - - /* The wavenumber is greater than the upper bound of the spectral range. - * Stop the copy since the remaining vertices are greater than the current - * one and are thus necessary out of the spectral range */ - if(nu >= spectral_range[1]) break; - - vtx->wavenumber = (float)nu; - vtx->ka = (float)ha; - ++i; - } - - /* Emit the last vertex for the upper bound of the spectral range if the line - * is clipped by it */ - if(line->wavenumber + mol_params->cutoff >= spectral_range[1]) { - struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices)+i; - const double nu = spectral_range[1]; - double ha = 0; - - switch(tree->args.mesh_type) { - case SLN_MESH_UPPER: - ha = nu > line->wavenumber - ? next_vertex_value(nu, wavenumbers, values) - : next_vertex_value(MIRROR(nu), wavenumbers, values); - break; - - case SLN_MESH_USUAL: - ha = line_value(tree, line, nu); - break; - - default: FATAL("Unreachable code.\n"); break; - } vtx->wavenumber = (float)nu; vtx->ka = (float)ha; - ++i; } #undef MIRROR - /* Several vertices could be skipped due to the clipping of the line by the - * spectral range. Resize the list of vertices to the effective number of - * registered vertices. */ - CHK(darray_vertex_resize(&tree->vertices, i) == RES_OK); + ASSERT(i == nvertices + line_nvertices); /* Setup the range of the line vertices */ vertices_range[0] = nvertices; - vertices_range[1] = i-1; + vertices_range[1] = i-1; /* Make the bound inclusive */ exit: return res; @@ -604,7 +501,7 @@ line_mesh (struct sln_tree* tree, const size_t iline, const size_t nvertices_hint, - size_t vertices_range[2]) /* out */ + size_t vertices_range[2]) /* out. Bounds are inclusive */ { /* The line */ struct line line = LINE_NULL; @@ -624,7 +521,7 @@ line_mesh darray_double_init(tree->sln->allocator, &wavenumbers); /* Setup the line wrt molecule concentration, isotope abundance, temperature - * and pression */ + * and pression */ res = line_setup(tree, iline, &line); if(res != RES_OK) goto error; diff --git a/src/sln_line.h b/src/sln_line.h @@ -66,7 +66,7 @@ line_mesh (struct sln_tree* tree, const size_t iline, const size_t nvertices_hint, - /* Range of generated line vertices. Upper bound is inclusive */ + /* Range of generated line vertices. Bounds are inclusive */ size_t vertices_range[2]); /* out */ #endif /* SLN_LINE_H */ diff --git a/src/sln_tree.c b/src/sln_tree.c @@ -351,7 +351,6 @@ sln_tree_create_from_stream READ(&tree->args.nvertices_hint, 1); READ(&tree->args.mesh_decimation_err, 1); READ(&tree->args.mesh_type, 1); - READ(&tree->args.wavenumber_range, 1); #undef READ res = shtr_isotope_metadata_create_from_stream(shtr, stream, &tree->args.metadata); @@ -465,7 +464,7 @@ sln_node_eval { double ka = 0; size_t iline; - ASSERT(tree && node && nu >= 0); + ASSERT(tree && node); FOR_EACH(iline, node->range[0], node->range[1]+1) { struct line line = LINE_NULL; @@ -490,7 +489,7 @@ sln_mesh_eval(const struct sln_mesh* mesh, const double wavenumber) const float nu = (float)wavenumber; size_t n; /* #vertices */ float u; /* Linear interpolation paraeter */ - ASSERT(mesh && nu >= 0 && mesh->nvertices); + ASSERT(mesh && mesh->nvertices); n = mesh->nvertices; @@ -547,7 +546,6 @@ sln_tree_write(const struct sln_tree* tree, FILE* stream) WRITE(&tree->args.nvertices_hint, 1); WRITE(&tree->args.mesh_decimation_err, 1); WRITE(&tree->args.mesh_type, 1); - WRITE(&tree->args.wavenumber_range, 1); #undef WRITE res = shtr_isotope_metadata_write(tree->args.metadata, stream); diff --git a/src/sln_tree_build.c b/src/sln_tree_build.c @@ -83,7 +83,7 @@ eval_ka { struct sln_mesh mesh = SLN_MESH_NULL; double ka = 0; - ASSERT(tree && node && wavenumber >= 0); + ASSERT(tree && node); switch(tree->args.mesh_type) { case SLN_MESH_UPPER: diff --git a/src/test_sln_tree.c b/src/test_sln_tree.c @@ -393,8 +393,6 @@ main(int argc, char** argv) tree_args.molecules[SHTR_CO2].cutoff = 50; tree_args.molecules[SHTR_O3 ].concentration = 1.0/3.0; tree_args.molecules[SHTR_O3 ].cutoff = 25; - tree_args.wavenumber_range[0] = 0; - tree_args.wavenumber_range[1] = INF; tree_args.pressure = 1; tree_args.temperature = 296;