star-line

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

commit 80ccb818abcc9d71d43deb824176265d9f44dcce
parent c414bb160ae6e65765a7dacbead3fd1fc6b21b1d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri,  6 May 2022 09:44:25 +0200

Add the nvertices_hit input argument to sln_tree_create

Internal refactoring of the line mesh function

Diffstat:
Msrc/sln.h | 7+++++--
Msrc/sln_line.c | 189++++++++++++++++++++++++++++++++++++++++++++++++++-----------------------------
Msrc/sln_line.h | 7+++++--
Msrc/sln_tree.c | 7+++++++
Msrc/sln_tree_build.c | 3+--
5 files changed, 138 insertions(+), 75 deletions(-)

diff --git a/src/sln.h b/src/sln.h @@ -105,11 +105,14 @@ struct sln_tree_create_args { size_t max_nlines_per_leaf; + /* Hint on the number of vertices around the the line center */ + size_t nvertices_hint; + enum sln_line_profile line_profile; }; #define SLN_TREE_CREATE_ARGS_DEFAULT__ \ - {NULL, NULL, {SLN_MOLECULE_NULL__}, 0, {0, DBL_MAX}, 0, 0, 64, \ - SLN_LINE_PROFILE_VOIGT} + {NULL, NULL, {SLN_MOLECULE_NULL__}, 0, {0, DBL_MAX}, 0, 0, 1, \ + 16, SLN_LINE_PROFILE_VOIGT} 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 @@ -32,8 +32,8 @@ #define AVOGADRO_NUMBER 6.02214076e23 /* molec.mol^-1 */ #define PERFECT_GAZ_CONSTANT 8.2057e-5 /* m^3.atm.mol^-1.K^-1 */ -#define MIN_NVERTICES_HINT 32 -#define MAX_NVERTICES_HINT 2048 +#define MIN_NVERTICES_HINT 8 +#define MAX_NVERTICES_HINT 128 STATIC_ASSERT(IS_POW2(MIN_NVERTICES_HINT), MIN_NVERTICES_HINT_is_not_a_pow2); STATIC_ASSERT(IS_POW2(MIN_NVERTICES_HINT), MAX_NVERTICES_HINT_is_not_a_pow2); @@ -178,7 +178,7 @@ error: * only the set of wavenumbers from the line center to its upper bound. */ static res_T regular_mesh_fragmented - (struct sln_tree* tree, + (const struct sln_tree* tree, const struct line* line, const struct molecule_params* mol_params, const double fragment_length, @@ -233,7 +233,7 @@ error: static res_T eval_mesh - (struct sln_tree* tree, + (const struct sln_tree* tree, const size_t iline, const struct darray_double* wavenumbers, struct darray_double* values) @@ -280,6 +280,112 @@ error: goto exit; } +/* Append the line mesh into the vertices array */ +static res_T +save_line_mesh + (const struct sln_tree* tree, + const size_t iline, + const struct molecule_params* mol_params, + const double spectral_range[2], + const struct darray_double* wavenumbers, + const struct darray_double* values, + struct darray_vertex* vertices) +{ + const struct line* line = NULL; + size_t nvertices; + size_t nwavenumbers; + size_t line_nvertices; + size_t ivertex; + size_t i; + res_T res = RES_OK; + ASSERT(tree && mol_params && wavenumbers && values && vertices); + ASSERT(darray_double_size_get(wavenumbers) == darray_double_size_get(values)); + ASSERT(spectral_range && spectral_range[0] <= spectral_range[1]); + + line = darray_line_cdata_get(&tree->lines) + iline; + nvertices = darray_vertex_size_get(vertices); + nwavenumbers = darray_double_size_get(wavenumbers); + + /* Compute the overall number of vertices of the line */ + line_nvertices = nwavenumbers + * 2 /* The line is symmetrical in its center */ + - 1;/* Do not duplicate the line center */ + + /* Allocate the list of line vertices */ + res = darray_vertex_resize(vertices, nvertices + line_nvertices); + if(res != RES_OK) { + log_err(tree->sln, "Error allocating the line vertices -- %s.\n", + res_to_cstr(res)); + goto error; + } + + i = nvertices; + + /* 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(vertices) + i; + vtx->wavenumber = (float)spectral_range[0]; + vtx->ka = (float)line_value(tree, iline, vtx->wavenumber); + ++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(vertices) + i; + const double ka = darray_double_cdata_get(values)[ivertex]; + const double nu = /* Mirror the wavenumber */ + 2*line->wavenumber - darray_double_cdata_get(wavenumbers)[ivertex]; + + /* The wavenumber is less than the lower bound of the spectral range. + * Discard the vertex */ + if(nu <= spectral_range[0]) + continue; + + vtx->wavenumber = (float)nu; + vtx->ka = (float)ka; + ++i; + } + + /* Copy the vertices of the line for its upper half */ + FOR_EACH(ivertex, 0, nwavenumbers) { + struct sln_vertex* vtx = darray_vertex_data_get(vertices) + i; + const double ka = darray_double_cdata_get(values)[ivertex]; + const double nu = darray_double_cdata_get(wavenumbers)[ivertex]; + + /* 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)ka; + ++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(vertices) + i; + vtx->wavenumber = (float)spectral_range[1]; + vtx->ka = (float)line_value(tree, iline, vtx->wavenumber); + ++i; + } + + /* 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(vertices, i) == RES_OK); + +exit: + return res; +error: + darray_vertex_resize(vertices, nvertices); + goto exit; +} + + /******************************************************************************* * Local function ******************************************************************************/ @@ -340,9 +446,10 @@ line_value res_T line_mesh - (struct sln_tree* tree, + (const struct sln_tree* tree, const size_t iline, - const size_t nvertices_hint) + const size_t nvertices_hint, + struct darray_vertex* vertices) { struct darray_double values; /* List of evaluated values */ struct darray_double wavenumbers; /* List of considered wavenumbers */ @@ -350,13 +457,8 @@ line_mesh const struct molecule_params* mol_params = NULL; const struct line* line = NULL; size_t nvertices_adjusted = 0; /* computed from nvertices_hint */ - size_t line_nvertices = 0; /* number of generated line vertices */ - size_t tree_nvertices = 0; /* #registered vertices of the tree */ - size_t ivertex; - size_t i; - res_T res = RES_OK; - ASSERT(tree && nvertices_hint); + ASSERT(tree && nvertices_hint && vertices); ASSERT(iline < darray_line_size_get(&tree->lines)); line = darray_line_cdata_get(&tree->lines) + iline; @@ -372,77 +474,26 @@ line_mesh (nvertices_hint, MIN_NVERTICES_HINT, MAX_NVERTICES_HINT); nvertices_adjusted = round_up_pow2(nvertices_adjusted); - /* Fetch the number already issued vertices */ - tree_nvertices = darray_vertex_size_get(&tree->vertices); - /* Retrieve the molecular parameters of the line to be mesh */ mol_params = tree_get_molecule_params(tree, shtr_line->molecule_id); + /* Emit the vertex coordinates, i.e. the wavenumbers */ res = regular_mesh_fragmented (tree, line, mol_params, line->gamma_l, nvertices_adjusted, &wavenumbers); if(res != RES_OK) goto error; + /* Evaluate the mesh vertices, i.e. define the line value for the list of + * wavenumbers */ eval_mesh(tree, iline, &wavenumbers, &values); - line_nvertices = darray_double_size_get(&wavenumbers); - res = darray_vertex_resize - (&tree->vertices, tree_nvertices + line_nvertices*2 - 1); - if(res != RES_OK) { - log_err(tree->sln, "Error allocating the line vertices -- %s.\n", - res_to_cstr(res)); - goto error; - } - - i = tree_nvertices; - - /* Emit the 1st vertex for the lower bound of the spectral range if the line - * cut it */ - if(line->wavenumber - mol_params->cutoff <= tree->wavenumber_range[0]) { - struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices) + i; - vtx->wavenumber = (float)tree->wavenumber_range[0]; - vtx->ka = (float)line_value(tree, iline, vtx->wavenumber); - ++i; - } - - /* Copy the first half of the line vertices */ - FOR_EACH_REVERSE(ivertex, line_nvertices-1, 0) { - struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices) + i; - double nu = (float)darray_double_cdata_get(&wavenumbers)[ivertex]; - nu = line->wavenumber - nu + line->wavenumber; - if(nu <= tree->wavenumber_range[0]) continue;/* Out of range */ - vtx->wavenumber = (float)nu; - vtx->ka = (float)darray_double_cdata_get(&values)[ivertex]; - ++i; - } - /* Copy the vertex of the line center and the second half of the line - * vertices */ - FOR_EACH(ivertex, 0, line_nvertices) { - struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices) + i; - const double nu = darray_double_cdata_get(&wavenumbers)[ivertex]; - if(nu >= tree->wavenumber_range[1]) break;/* Out of range */ - vtx->wavenumber = (float)nu; - vtx->ka = (float)darray_double_cdata_get(&values)[ivertex]; - ++i; - } - - /* Emit the last vertex for the upper bound of the spectral range */ - if(line->wavenumber + mol_params->cutoff >= tree->wavenumber_range[1]) { - struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices) + i; - vtx->wavenumber = (float)tree->wavenumber_range[1]; - vtx->ka = (float)line_value(tree, iline, vtx->wavenumber); - ++i; - } - - /* Several vertices could be deleted due to spectral range. So, resize the - * list of vertices to the effective number of registered vertices. */ - CHK(darray_vertex_resize(&tree->vertices, i) == RES_OK); + res = save_line_mesh(tree, iline, mol_params, tree->wavenumber_range, + &wavenumbers, &values, vertices); + if(res != RES_OK) goto error; exit: darray_double_release(&values); darray_double_release(&wavenumbers); return res; error: - /* Clean up registered vertices */ - darray_vertex_resize(&tree->vertices, tree_nvertices); goto exit; } diff --git a/src/sln_line.h b/src/sln_line.h @@ -32,6 +32,7 @@ struct line { /* Forward declaration */ struct sln_tree; +struct darray_vertex; static INLINE double line_center @@ -53,10 +54,12 @@ line_value const size_t iline, const double wavenumber); +/* Append the line vertices into `vertices' */ extern LOCAL_SYM res_T line_mesh - (struct sln_tree* tree, + (const struct sln_tree* tree, const size_t iline, - const size_t nvertices_hint); + const size_t nvertices_hint, + struct darray_vertex* vertices); #endif /* SLN_LINE_H */ diff --git a/src/sln_tree.c b/src/sln_tree.c @@ -194,6 +194,13 @@ check_sln_tree_create_args return RES_BAD_ARG; } + if(args->nvertices_hint == 0) { + log_err(sln, + "%s: invalid hint on the number of vertices around the line center %lu.\n", + func_name, (unsigned long)args->nvertices_hint); + return RES_BAD_ARG; + } + if((unsigned)args->line_profile >= SLN_LINE_PROFILES_COUNT__) { log_err(sln, "%s: invalid line profile %d.\n", func_name, args->line_profile); diff --git a/src/sln_tree_build.c b/src/sln_tree_build.c @@ -40,7 +40,6 @@ setup_leaf struct sln_node* leaf) { res_T res = RES_OK; - (void)args; /* It should define the number of vertices */ /* Currently assume that we have only one line per leaf */ ASSERT(leaf->range[0] == leaf->range[1]); @@ -49,7 +48,7 @@ setup_leaf leaf->ivertices = darray_vertex_size_get(&tree->vertices); - res = line_mesh(tree, leaf->range[0], 64/*TODO it must be an argument*/); + res = line_mesh(tree, leaf->range[0], args->nvertices_hint, &tree->vertices); if(res != RES_OK) goto error; leaf->nvertices = ui64_to_ui32