commit 7da5293205e492c5fb7eb96c200d29ab7ac21d7d
parent 144a2ac7d11d8a9a31021d36450577ce75c19006
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 22 Sep 2026 11:46:04 +0200
Let the caller change the thermo properties of a tree's lines
The sln_node_eval and sln_node_get_line functions now take an optional
structured argument that defines the temperature, pressure and the
species concentration of the node line(s). If the corresponding argument
argument is NULL it is the thermodynamic poprieties to which the tree
was built that are then used.
This API modification allows the caller to use the same tree to sample
lines at different thermodynmic conditions than those used to build the
tree. The user must nevertheless ensure that the criteria for sampling
by importance are met, in particular that the polyline of a node is
an upper bound of the sampled line(s).
Implementation is careful not to increase the cost of the functions that
this commit updates. This could be the case if the concentration list was
checked at each function call. To avoid this extra cost, the
sln_node_get_line function only checks the concentration of the single
molecule corresponding to the input line index. While sln_node_eval
checks the consistency of concentrations but only in debug.
No tests were performed. Only the default behavior was verified to give
the same results as before.
Diffstat:
8 files changed, 219 insertions(+), 63 deletions(-)
diff --git a/src/sln.h b/src/sln.h
@@ -210,6 +210,15 @@ struct sln_tree_desc {
}
static const struct sln_tree_desc SLN_TREE_DESC_NULL = SLN_TREE_DESC_NULL__;
+struct sln_thermo_props {
+ double concentrations[SHTR_MAX_MOLECULE_COUNT];
+ double pressure; /* [atm] */
+ double temperature; /* [K] */
+};
+#define SLN_THERMO_PROPS_NULL__ {{0},0,0}
+static const struct sln_thermo_props SLN_THERMO_PROPS_NULL =
+ SLN_THERMO_PROPS_NULL__;
+
struct sln_node_desc {
/* Range of lines belonging to the node. The endpoints are included */
size_t ilines[2];
@@ -351,6 +360,9 @@ SLN_API res_T
sln_tree_get_line
(const struct sln_tree* tree,
const size_t iline,
+ /* Thermodynamic properties to which the line is recovered.
+ * Can be NULL, so these properties are those used to build the tree */
+ const struct sln_thermo_props* props,
struct sln_line* line);
SLN_API res_T
@@ -381,6 +393,9 @@ SLN_API double
sln_node_eval
(const struct sln_tree* tree,
const struct sln_node* node,
+ /* Thermodynamic properties to which the node lines are evaluated.
+ * Can be NULL, so these properties are those used to build the tree */
+ const struct sln_thermo_props* props,
const double wavenumber); /* In cm^-1 */
SLN_API res_T
diff --git a/src/sln_get.c b/src/sln_get.c
@@ -427,7 +427,7 @@ print_node_value(const struct cmd* cmd)
if(res != RES_OK) goto error;
val_mesh = sln_mesh_eval(&mesh, cmd->args.wavenumber);
- val_node = sln_node_eval(cmd->tree, node, cmd->args.wavenumber);
+ val_node = sln_node_eval(cmd->tree, node, NULL, cmd->args.wavenumber);
printf("ka(%e) = %e ~ %e\n", cmd->args.wavenumber, val_node, val_mesh);
diff --git a/src/sln_line.c b/src/sln_line.c
@@ -427,12 +427,15 @@ res_T
line_setup
(const struct sln_tree* tree,
const size_t iline,
+ const struct sln_thermo_props* props,
struct sln_line* line)
{
struct shtr_molecule molecule = SHTR_MOLECULE_NULL;
struct shtr_line shtr_line = SHTR_LINE_NULL;
- double molar_mass = 0; /* In kg.mol^-1 */
- const struct sln_molecule* mol_params = NULL;
+ double concentration = 0;
+ double molar_mass = 0; /*[kg.mol^-1]*/
+ double pressure = 0; /*[atm]*/
+ double temperature = 0; /*[K]*/
res_T res = RES_OK;
ASSERT(tree && line);
@@ -443,7 +446,15 @@ line_setup
ASSERT(!SHTR_MOLECULE_IS_NULL(&molecule));
ASSERT(molecule.nisotopes > (size_t)shtr_line.isotope_id_local);
- mol_params = tree->args.molecules + shtr_line.molecule_id;
+ if(!props) { /* Use thermo properties used to build the tree */
+ concentration = tree->args.molecules[shtr_line.molecule_id].concentration;
+ pressure = tree->args.pressure; /*[atm]*/
+ temperature = tree->args.temperature; /*[K]*/
+ } else {
+ concentration = props->concentrations[shtr_line.molecule_id];
+ pressure = props->pressure; /*[atm]*/
+ temperature = props->temperature; /*[K]*/
+ }
/* Convert the molar mass of the line from g.mol^-1 to kg.mol^-1 */
molar_mass = molecule.isotopes[shtr_line.isotope_id_local].molar_mass*1e-3;
@@ -452,12 +463,12 @@ line_setup
res = line_profile_factor(tree, &shtr_line, &line->profile_factor);
if(res != RES_OK) goto error;
- line->wavenumber = line_center(&shtr_line, tree->args.pressure);
+ line->wavenumber = line_center(&shtr_line, pressure);
line->gamma_d = sln_compute_line_half_width_doppler
- (line->wavenumber, molar_mass, tree->args.temperature);
+ (line->wavenumber, molar_mass, temperature);
line->gamma_l = sln_compute_line_half_width_lorentz
- (shtr_line.gamma_air, shtr_line.gamma_self, tree->args.temperature,
- tree->args.pressure, shtr_line.n_air, mol_params->concentration);
+ (shtr_line.gamma_air, shtr_line.gamma_self, temperature,
+ pressure, shtr_line.n_air, concentration);
line->molecule_id = shtr_line.molecule_id;
exit:
@@ -493,7 +504,7 @@ line_mesh
/* Setup the line wrt molecule concentration, isotope abundance, temperature
* and pressure */
- res = line_setup(tree, iline, &line);
+ res = line_setup(tree, iline, NULL/*default thermo props*/, &line);
if(res != RES_OK) goto error;
/* Adjust the hint on the number of vertices. This is not actually the real
diff --git a/src/sln_line.h b/src/sln_line.h
@@ -43,6 +43,9 @@ extern LOCAL_SYM res_T
line_setup
(const struct sln_tree* tree,
const size_t iline,
+ /* Thermodynamic properties to which the line is recovered.
+ * Can be NULL, so these properties are those used to build the tree */
+ const struct sln_thermo_props* props,
struct sln_line* line);
extern LOCAL_SYM res_T
diff --git a/src/sln_slab.c b/src/sln_slab.c
@@ -411,7 +411,7 @@ realisation
/* Evaluate the value of the line and compute the probability of being
* absorbed by it */
- leaf_ka = sln_node_eval(cmd->tree, leaf, nu);
+ leaf_ka = sln_node_eval(cmd->tree, leaf, NULL, nu);
proba_abs = leaf_ka / (leaf_proba*ka_max);
if((res = check_proba(cmd, proba_abs)) != RES_OK) goto error;
diff --git a/src/sln_tree.c b/src/sln_tree.c
@@ -40,20 +40,52 @@ static const struct stream STREAM_NULL = {NULL, NULL, 0};
/*******************************************************************************
* Helper functions
******************************************************************************/
-/* Check that the sum of the molecular concentrations is equal to 1 */
-static res_T
+static INLINE res_T
check_molecule_concentration
- (struct sln_device* sln,
+ (const struct sln_device* sln,
const char* caller,
- const struct sln_tree_create_args* args)
+ const enum shtr_molecule_id molecule_id,
+ const double concentration)
+{
+ ASSERT(sln && caller);
+
+ if(concentration == 0) {
+ /* A molecular concentration of zero is allowed, but may be a user error,
+ * as 0 is the default concentration in the tree creation arguments.
+ * Therefore, warn the user about this value so that they can determine
+ * whether or not it is an error on their part. */
+ WARN(sln, "%s: the concentration of %s is zero.\n",
+ caller, shtr_molecule_cstr(molecule_id));
+
+ } else if(concentration < 0) {
+ /* Concentration cannot be negative... */
+ ERROR(sln, "%s: invalid %s concentration: %g.\n",
+ FUNC_NAME, shtr_molecule_cstr(molecule_id),
+ concentration);
+ return RES_BAD_ARG;
+ }
+
+ return RES_OK;
+}
+
+static res_T
+check_concentrations_list
+ (const struct sln_device* sln,
+ const char* caller,
+ const double concentrations[SHTR_MAX_MOLECULE_COUNT])
{
- int i = 0;
double sum = 0;
- ASSERT(sln && caller && args);
+ int i = 0;
+ res_T res = RES_OK;
+ ASSERT(sln && caller && concentrations);
FOR_EACH(i, 0, SHTR_MAX_MOLECULE_COUNT) {
if(i == SHTR_MOLECULE_ID_NULL) continue;
- sum += args->molecules[i].concentration;
+
+ res = check_molecule_concentration(sln, caller, i, concentrations[i]);
+ if(res != RES_OK) goto error;
+
+ sum += concentrations[i];
}
/* The sum of molecular concentrations must be less than or equal to 1. It may
@@ -63,10 +95,32 @@ check_molecule_concentration
ERROR(sln,
"%s: the sum of molecule concentrations is greater than 1: %g\n",
caller, sum);
- return RES_BAD_ARG;
+ res = RES_BAD_ARG;
+ goto error;
}
- return RES_OK;
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/* Check the consistency of the molecular concentrations */
+static INLINE res_T
+check_mixture_concentrations
+ (struct sln_device* sln,
+ const char* caller,
+ const struct sln_tree_create_args* args)
+{
+ double concentrations[SHTR_MAX_MOLECULE_COUNT] = {0};
+ int i = 0;
+ ASSERT(sln && caller && args);
+
+ FOR_EACH(i, 0, SHTR_MAX_MOLECULE_COUNT) {
+ concentrations[i] = args->molecules[i].concentration;
+ }
+
+ return check_concentrations_list(sln, caller, concentrations);
}
/* Verify that the isotope abundance are valids */
@@ -115,14 +169,13 @@ check_molecules
{
char molecule_ok[SHTR_MAX_MOLECULE_COUNT] = {0};
- double concentrations_sum = 0;
size_t iline = 0;
size_t nlines = 0;
res_T res = RES_OK;
ASSERT(args->lines);
- res = check_molecule_concentration(sln, caller, args);
- if(res != RES_OK) return res;
+ res = check_mixture_concentrations(sln, caller, args);
+ if(res != RES_OK) goto error;
/* Iterate over the lines to define which molecules has to be checked, i.e.,
* the ones used in the mixture */
@@ -138,24 +191,6 @@ check_molecules
molecule = args->molecules + line.molecule_id;
- if(molecule->concentration == 0) {
- /* A molecular concentration of zero is allowed, but may be a user error,
- * as 0 is the default concentration in the tree creation arguments.
- * Therefore, warn the user about this value so that they can determine
- * whether or not it is an error on their part. */
- WARN(sln, "%s: the concentration of %s is zero.\n",
- caller, shtr_molecule_cstr(line.molecule_id));
-
- } else if(molecule->concentration < 0) {
- /* Concentration cannot be negative... */
- ERROR(sln, "%s: invalid %s concentration: %g.\n",
- FUNC_NAME, shtr_molecule_cstr(line.molecule_id),
- molecule->concentration);
- return RES_BAD_ARG;
- }
-
- concentrations_sum += molecule->concentration;
-
if(molecule->cutoff <= 0) {
/* ... cutoff either */
ERROR(sln, "%s: invalid %s cutoff: %g.\n",
@@ -164,21 +199,40 @@ check_molecules
}
res = check_molecule_isotope_abundances(sln, caller, molecule);
- if(res != RES_OK) return res;
+ if(res != RES_OK) goto error;
molecule_ok[line.molecule_id] = 1;
}
- /* The sum of molecular concentrations must be less than or equal to 1. It may
- * be less than 1 if the remaining part of the mixture is (implicitly) defined
- * as a radiatively inactive gas */
- if(concentrations_sum > 1 && (concentrations_sum - 1) > 1e-6) {
- ERROR(sln,
- "%s: the sum of molecule concentrations is greater than 1: %g\n",
- caller, concentrations_sum);
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static INLINE res_T
+check_pressure
+ (const struct sln_device* sln,
+ const char* caller,
+ const double pressure /*[atm]*/)
+{
+ if(pressure < 0) {
+ ERROR(sln, "%s: invalid negative pressure %g atm\n", caller, pressure);
return RES_BAD_ARG;
}
+ return RES_OK;
+}
+static INLINE res_T
+check_temperature
+ (const struct sln_device* sln,
+ const char* caller,
+ const double temperature /*[K]*/)
+{
+ if(temperature < 0) {
+ ERROR(sln, "%s: invalid negative temperature %g K\n", caller, temperature);
+ return RES_BAD_ARG;
+ }
return RES_OK;
}
@@ -203,6 +257,14 @@ check_sln_tree_create_args
return RES_BAD_ARG;
}
+ if((res = check_pressure(sln, caller, args->pressure)) != RES_OK) {
+ return res;
+ }
+
+ if((res = check_temperature(sln, caller, args->temperature)) != RES_OK) {
+ return res;
+ }
+
if(args->nvertices_hint == 0) {
ERROR(sln,
"%s: invalid hint on the number of vertices around the line center %lu.\n",
@@ -296,6 +358,64 @@ check_sln_tree_write_args
return RES_OK;
}
+
+static res_T
+check_line_thermo_props
+ (const struct sln_tree* tree,
+ const char* caller,
+ const size_t iline,
+ const struct sln_thermo_props* props)
+{
+ struct shtr_line line = SHTR_LINE_NULL;
+ res_T res = RES_OK;
+ ASSERT(tree && caller);
+
+ if(!props) goto exit; /* Default thermo props */
+
+ SHTR(line_list_at(tree->args.lines, iline, &line));
+
+ res = check_molecule_concentration(tree->sln, caller, line.molecule_id,
+ props->concentrations[line.molecule_id]);
+ if(res != RES_OK) goto error;
+
+ res = check_pressure(tree->sln, caller, props->pressure);
+ if(res != RES_OK) goto error;
+
+ res = check_temperature(tree->sln, caller, props->temperature);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static INLINE res_T
+check_sln_thermo_props
+ (const struct sln_device* sln,
+ const char* caller,
+ const struct sln_thermo_props* props)
+{
+ res_T res = RES_OK;
+ ASSERT(sln && caller);
+
+ if(!props) goto exit; /* Default thermo props */
+
+ res = check_concentrations_list(sln, caller, props->concentrations);
+ if(res != RES_OK) goto error;
+
+ res = check_pressure(sln, caller, props->pressure);
+ if(res != RES_OK) goto error;
+
+ res = check_temperature(sln, caller, props->temperature);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
static INLINE void
stream_release(struct stream* stream)
{
@@ -703,6 +823,7 @@ res_T
sln_tree_get_line
(const struct sln_tree* tree,
const size_t iline,
+ const struct sln_thermo_props* props,
struct sln_line* line)
{
size_t nlines = 0;
@@ -713,7 +834,10 @@ sln_tree_get_line
SHTR(line_list_get_size(tree->args.lines, &nlines));
if(iline >= nlines) { res = RES_BAD_ARG; goto error; }
- res = line_setup(tree, iline, line);
+ res = check_line_thermo_props(tree, FUNC_NAME, iline, props);
+ if(res != RES_OK) goto error;
+
+ res = line_setup(tree, iline, props, line);
if(res != RES_OK) {
ERROR(tree->sln, "%s: could not setup the line %lu-- %s\n",
FUNC_NAME, iline, res_to_cstr(res));
@@ -775,17 +899,19 @@ double
sln_node_eval
(const struct sln_tree* tree,
const struct sln_node* node,
+ const struct sln_thermo_props* props,
const double nu)
{
double ka = 0;
size_t iline;
ASSERT(tree && node);
+ ASSERT(check_sln_thermo_props(tree->sln, FUNC_NAME, props) == RES_OK);
FOR_EACH(iline, node->range[0], node->range[1]+1) {
struct sln_line line = SLN_LINE_NULL;
res_T res = RES_OK;
- res = line_setup(tree, iline, &line);
+ res = line_setup(tree, iline, props, &line);
if(res != RES_OK) {
WARN(tree->sln, "%s: could not setup the line %lu-- %s\n",
FUNC_NAME, iline, res_to_cstr(res));
diff --git a/src/sln_tree_build.c b/src/sln_tree_build.c
@@ -270,8 +270,9 @@ build_leaf_polyline_from_Nlines
const size_t iline = leaf->range[0] + i;
/* Mesh the line in the temporary vertex buffer */
- res = line_mesh(tree, iline, tree->args.nvertices_hint, &scratch->vertices,
- vertices_range);
+ res = line_mesh
+ (/* in */ tree, iline, tree->args.nvertices_hint,
+ /* out */ &scratch->vertices, vertices_range);
if(res != RES_OK) goto error;
/* Decimate the line mesh */
diff --git a/src/test_sln_tree.c b/src/test_sln_tree.c
@@ -179,8 +179,8 @@ check_node_equality
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(sln_tree_get_line(tree1, iline, NULL, &line1) == RES_OK);
+ CHK(sln_tree_get_line(tree2, iline, NULL, &line2) == RES_OK);
CHK(line1.wavenumber == line2.wavenumber);
CHK(line1.profile_factor == line2.profile_factor);
@@ -267,15 +267,15 @@ check_node_value(const struct sln_tree* tree, const struct sln_node* node)
CHK(sln_node_get_desc(tree, node, &desc) == RES_OK);
- CHK(sln_tree_get_line(tree, desc.ilines[0], &line) == RES_OK);
+ CHK(sln_tree_get_line(tree, desc.ilines[0], NULL, &line) == RES_OK);
nu = line.wavenumber + 10;
- CHK(sln_tree_get_line(tree, desc.ilines[1], &line) == RES_OK);
+ CHK(sln_tree_get_line(tree, desc.ilines[1], NULL, &line) == RES_OK);
nu += line.wavenumber - 10;
nu *= 0.5;
- ka_node = sln_node_eval(tree, node, nu);
+ ka_node = sln_node_eval(tree, node, NULL, nu);
FOR_EACH(iline, desc.ilines[0], desc.ilines[1]+1/*inclusive*/) {
- CHK(sln_tree_get_line(tree, iline, &line) == RES_OK);
+ CHK(sln_tree_get_line(tree, iline, NULL, &line) == RES_OK);
ka_ref += sln_line_eval(tree, &line, nu);
}
@@ -329,11 +329,11 @@ 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(sln_tree_get_line(NULL, 0, NULL, &line) == RES_BAD_ARG);
+ CHK(sln_tree_get_line(tree, nlines, NULL, &line) == RES_BAD_ARG);
+ CHK(sln_tree_get_line(tree, 0, NULL, NULL) == RES_BAD_ARG);
+ CHK(sln_tree_get_line(tree, 0, NULL, &line) == RES_OK);
+ CHK(sln_tree_get_line(tree, nlines-1, NULL, &line) == RES_OK);
CHK(node = sln_tree_get_root(tree));
CHK(node != NULL);
