commit a9389fa526a6e80edb6e4953d00deeb0e6c52ce1
parent 8f0110bb96b46de17409c62e43bb4bd00c172048
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 3 May 2022 16:28:51 +0200
Precompute per line data
This will speed up the evaluation of the line value.
Diffstat:
5 files changed, 474 insertions(+), 11 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -29,12 +29,16 @@ option(NO_TEST "Do not build tests" OFF)
find_package(RCMake 0.4 REQUIRED)
find_package(RSys 0.12.1 REQUIRED)
find_package(StarHITRAN REQUIRED)
+find_package(LBLU REQUIRED)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
include(rcmake)
include(rcmake_runtime)
-include_directories(${RSys_INCLUDE_DIR} ${StarHITRAN_INCLUDE_DIR})
+include_directories(
+ ${LBLU_INCLUDE_DIR}
+ ${RSys_INCLUDE_DIR}
+ ${StarHITRAN_INCLUDE_DIR})
################################################################################
# Configure and define targets
@@ -46,12 +50,14 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SLN_FILES_SRC
sln_device.c
+ sln_line.c
sln_log.c
sln_tree.c
sln_tree_build.c
sln_voigt.c)
set(SLN_FILES_INC
sln_device_c.h
+ sln_line.h
sln_log.h)
set(SLN_FILES_INC_API
sln.h)
@@ -66,7 +72,7 @@ rcmake_prepend_path(SLN_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
add_library(sln SHARED
${SLN_FILES_SRC} ${SLN_FILES_INC} ${SLN_FILES_INC_API})
-target_link_libraries(sln RSys StarHITRAN m)
+target_link_libraries(sln LBLU RSys StarHITRAN m)
set_target_properties(sln PROPERTIES
DEFINE_SYMBOL SLN_SHARED_BUILD
diff --git a/src/sln_line.c b/src/sln_line.c
@@ -0,0 +1,337 @@
+/* Copyright (C) 2022 CNRS - LMD
+ * Copyright (C) 2022 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2022 Université Paul Sabatier - IRIT
+ * Copyright (C) 2022 Université Paul Sabatier - Laplace
+ *
+ * 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 "sln_device_c.h"
+#include "sln_line.h"
+#include "sln_log.h"
+#include "sln_tree_c.h"
+
+#include <lblu.h>
+#include <star/shtr.h>
+
+#include <rsys/cstr.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/math.h>
+
+#define T_REF 296.0 /* K */
+#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
+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);
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static INLINE double
+line_intensity
+ (const double intensity_ref, /* Reference intensity [cm^-1/(molec.cm^2)] */
+ const double lower_state_energy, /* [cm^-1] */
+ const double partition_function,
+ const double temperature, /* [K] */
+ const double temperature_ref, /* [K] */
+ const double wavenumber) /* [cm^-1] */
+{
+ const double C2 = 1.4388; /* 2nd Planck constant [K.cm] */
+
+ const double fol = /* TODO ask to Yaniss why this variable is named fol */
+ (1-exp(-C2*wavenumber/temperature))
+ / (1-exp(-C2*wavenumber/temperature_ref));
+
+ const double tmp =
+ exp(-C2*lower_state_energy/temperature)
+ / exp(-C2*lower_state_energy/temperature_ref);
+
+ return intensity_ref * partition_function * tmp * fol ;
+}
+
+static double
+line_profile_factor
+ (const struct sln_tree* tree,
+ const struct shtr_line* shtr_line)
+{
+ /* Star-HITRAN data */
+ struct shtr_molecule molecule = SHTR_MOLECULE_NULL;
+ const struct shtr_isotope* isotope = NULL;
+
+ /* Mixture parameters */
+ const struct molecule_params* mol_params = NULL;
+
+ /* Miscellaneous */
+ double density; /* In molec.cm^-3 */
+ double intensity, intensity_ref; /* In cm^-1/(molec.cm^2) */
+ double Q, Q_T, Q_Tref; /* Partition function */
+ double nu_c; /* In cm^-1 */
+ double profile_factor; /* In m^-1.cm^-1 */
+ double gj; /* State independant degeneracy factor */
+ double Ps; /* In atm */
+ double T; /* Temperature */
+ int molid; /* Molecule id */
+ int isoid; /* Isotope id local to its molecule */
+ ASSERT(tree);
+
+ /* Fetch the molecule data */
+ mol_params = tree_get_molecule_params(tree, shtr_line->molecule_id);
+ SHTR(isotope_metadata_find_molecule
+ (tree->metadata, shtr_line->molecule_id, &molecule));
+ ASSERT(!SHTR_MOLECULE_IS_NULL(&molecule));
+ ASSERT(molecule.nisotopes > (size_t)shtr_line->isotope_id_local);
+ isotope = molecule.isotopes + shtr_line->isotope_id_local;
+
+ nu_c = line_center(shtr_line, tree->pressure);
+
+ /* Compute the intensity */
+ Ps = tree->pressure * mol_params->concentration;
+ density = (AVOGADRO_NUMBER * Ps) / (PERFECT_GAZ_CONSTANT * tree->temperature);
+ density = density * 1e-6; /* Convert in molec.cm^-3 */
+
+ /* Compute the partition function TODO check that Q296K is equal to the
+ * partition function returned by BD_TIPS_2017(T_REF) */
+ Q_Tref = isotope->Q296K;
+ molid = shtr_line->molecule_id;
+ isoid = shtr_line->isotope_id_local;
+ T = tree->temperature;
+ BD_TIPS_2017(&molid, &T, &isoid, &gj, &Q_T);
+ Q = Q_Tref/Q_T;
+
+ /* Compute the intensity */
+ intensity_ref =
+ shtr_line->intensity/isotope->abundance
+ * mol_params->isotopes_abundance[shtr_line->isotope_id_local];
+ intensity = line_intensity(intensity_ref, shtr_line->lower_state_energy, Q,
+ tree->temperature, T_REF, nu_c);
+
+ profile_factor = 1.e2 * density * intensity; /* In m^-1.cm^-1 */
+ return profile_factor;
+}
+
+/* Regularly mesh the interval [wavenumber, wavenumber+spectral_length[. Note
+ * that the upper bound is excluded, this means that the last vertex of the
+ * interval is not emitted */
+static INLINE res_T
+regular_mesh
+ (const double wavenumber, /* Wavenumber where the mesh begins [cm^-1] */
+ const double spectral_length, /* Size of the spectral interval to mesh [cm^-1] */
+ const size_t nvertices, /* #vertices to issue */
+ struct darray_double* wavenumbers) /* List of issued vertices */
+{
+ /* Do not issue the vertex on the upper bound of the spectral range. That's
+ * why we assume that the number of steps is equal to the number of vertices
+ * and not to the number of vertices minus 1 */
+ const double step = spectral_length / (double)nvertices;
+ size_t ivtx;
+ res_T res = RES_OK;
+ ASSERT(wavenumber >= 0 && spectral_length > 0 && wavenumbers);
+
+ FOR_EACH(ivtx, 0, nvertices) {
+ const double nu = wavenumber + (double)ivtx*step;
+ res = darray_double_push_back(wavenumbers, &nu);
+ if(res != RES_OK) goto error;
+ }
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/* The line is regularly discretized into a set of fragments of fixed size.
+ * Their discretization is finer for the fragments around the center of the
+ * line and becomes coarser as the fragments move away from it. Note that a
+ * line is symmetrical in its center. As a consequence, the returned list is
+ * 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 line* line,
+ const struct molecule_params* mol_params,
+ const double fragment_length,
+ const size_t nvertices,
+ struct darray_double* wavenumbers) /* List of issued vertices */
+{
+ /* Fragment parameters */
+ double fragment_nu_min = 0; /* Lower bound of the fragment */
+ size_t fragment_nvtx = 0; /* #vertices into the fragment */
+
+ /* Miscellaneous */
+ double line_nu_min = 0; /* In cm^-1 */
+ double line_nu_max = 0; /* In cm^-1 */
+ res_T res = RES_OK;
+
+ ASSERT(tree && line && mol_params && wavenumbers);
+ ASSERT(fragment_length > 0);
+ ASSERT(IS_POW2(nvertices));
+
+ /* Compute the spectral range of the line from its center to its cutoff */
+ line_nu_min = line->wavenumber;
+ line_nu_max = line->wavenumber + mol_params->cutoff;
+
+ /* Define the number of vertices for the first interval in [nu, gamma_l] */
+ fragment_nu_min = line_nu_min;
+ fragment_nvtx = MMAX(nvertices/2, 2);
+
+ while(fragment_nu_min < line_nu_max) {
+ const double spectral_length =
+ MMIN(fragment_length, mol_params->cutoff - fragment_nu_min);
+
+ res = regular_mesh
+ (fragment_nu_min, spectral_length, fragment_nvtx, wavenumbers);
+ if(res != RES_OK) goto error;
+
+ fragment_nu_min += fragment_length;
+ fragment_nvtx = MMAX(fragment_nvtx/2, 2);
+ }
+
+ /* Register the last vertex, i.e. the upper bound of the spectral range */
+ res = darray_double_push_back(wavenumbers, &line_nu_max);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ log_err(tree->sln, "Error meshing the line -- %s.\n", res_to_cstr(res));
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+void
+line_setup(struct sln_tree* tree, const size_t iline)
+{
+ struct shtr_molecule molecule = SHTR_MOLECULE_NULL;
+ const struct shtr_line* shtr_line = NULL;
+ struct line* line = NULL;
+ double molar_mass = 0; /* In kg.mol^-1 */
+ const struct molecule_params* mol_params = NULL;
+ ASSERT(tree && iline < darray_line_size_get(&tree->lines));
+
+ line = darray_line_data_get(&tree->lines) + iline;
+
+ SHTR(lines_view_get_line(tree->lines_view, iline, &shtr_line));
+ SHTR(isotope_metadata_find_molecule
+ (tree->metadata, shtr_line->molecule_id, &molecule));
+ ASSERT(!SHTR_MOLECULE_IS_NULL(&molecule));
+ ASSERT(molecule.nisotopes > (size_t)shtr_line->isotope_id_local);
+ mol_params = tree_get_molecule_params(tree, shtr_line->molecule_id);
+
+ /* 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;
+
+ /* Setup the line */
+ line->wavenumber = line_center(shtr_line, tree->pressure);
+ line->profile_factor = line_profile_factor(tree, shtr_line);
+ line->gamma_d = sln_compute_line_half_width_doppler
+ (line->wavenumber, molar_mass, tree->temperature);
+ line->gamma_l = sln_compute_line_half_width_lorentz(shtr_line->gamma_air,
+ shtr_line->gamma_self, tree->pressure, mol_params->concentration);
+}
+
+double
+line_value
+ (const struct sln_tree* tree,
+ const size_t iline,
+ const double wavenumber)
+{
+ const struct line* line = NULL;
+ double profile;
+ ASSERT(tree && iline < darray_line_size_get(&tree->lines));
+
+ line = darray_line_cdata_get(&tree->lines) + iline;
+ profile = sln_compute_voigt_profile
+ (wavenumber, line->wavenumber, line->gamma_d, line->gamma_l);
+ return line->profile_factor * profile;
+}
+
+res_T
+line_mesh
+ (struct sln_tree* tree,
+ const size_t iline,
+ const size_t nvertices_hint)
+{
+ struct darray_double wavenumbers; /* List of considered wavenumbers */
+ const struct shtr_line* shtr_line = NULL;
+ 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 && line && nvertices_hint);
+ ASSERT(iline < darray_line_size_get(&tree->lines));
+
+ line = darray_line_cdata_get(&tree->lines) + iline;
+ SHTR(lines_view_get_line(tree->lines_view, iline, &shtr_line));
+
+ darray_double_init(tree->sln->allocator, &wavenumbers);
+
+ /* Adjust the hint on the number of vertices. This is not actually the real
+ * number of vertices but a adjusted hint on it. This new value ensures that
+ * it is a power of 2 included in [MIN_NVERTICES_HINT, MAX_NVERTICES_HINT] */
+ nvertices_adjusted = CLAMP
+ (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);
+
+ res = regular_mesh_fragmented
+ (tree, line, mol_params, line->gamma_l, nvertices_adjusted, &wavenumbers);
+ if(res != RES_OK) goto error;
+
+ 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;
+ }
+
+ /* TODO compute the line value for each vertex */
+
+ /* Copy the first half of the line vertices */
+ i = tree_nvertices;
+ FOR_EACH_REVERSE(ivertex, line_nvertices-1, 0) {
+ darray_vertex_data_get(&tree->vertices)[i].wavenumber =
+ (float)darray_double_cdata_get(&wavenumbers)[ivertex];
+ ++i;
+ }
+ /* Copy the vertex of the line center and the second half of the line
+ * vertices */
+ FOR_EACH(ivertex, 0, line_nvertices) {
+ darray_vertex_data_get(&tree->vertices)[i].wavenumber =
+ (float)darray_double_cdata_get(&wavenumbers)[ivertex];
+ ++i;
+ }
+
+exit:
+ 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
@@ -0,0 +1,62 @@
+/* Copyright (C) 2022 CNRS - LMD
+ * Copyright (C) 2022 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2022 Université Paul Sabatier - IRIT
+ * Copyright (C) 2022 Université Paul Sabatier - Laplace
+ *
+ * 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/>. */
+
+#ifndef SLN_LINE_H
+#define SLN_LINE_H
+
+#include <star/shtr.h>
+#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 */
+};
+
+/* Forward declaration */
+struct sln_tree;
+
+static INLINE double
+line_center
+ (const struct shtr_line* line,
+ const double pressure) /* In atm */
+{
+ ASSERT(line && pressure >= 0);
+ return line->wavenumber + line->delta_air * pressure;
+}
+
+extern LOCAL_SYM void
+line_setup
+ (struct sln_tree* tree,
+ const size_t iline);
+
+extern LOCAL_SYM double
+line_value
+ (const struct sln_tree* tree,
+ const size_t iline,
+ const double wavenumber);
+
+extern LOCAL_SYM res_T
+line_mesh
+ (struct sln_tree* tree,
+ const size_t iline,
+ const size_t nvertices_hint);
+
+#endif /* SLN_LINE_H */
diff --git a/src/sln_tree.c b/src/sln_tree.c
@@ -302,9 +302,41 @@ store_input_args
tree->wavenumber_range[0] = args->wavenumber_range[0];
tree->wavenumber_range[1] = args->wavenumber_range[1];
tree->max_nlines_per_leaf = args->max_nlines_per_leaf;
+ tree->temperature = args->temperature;
+ tree->pressure = args->pressure;
+ SHTR(isotope_metadata_ref_get(args->metadata));
+ tree->metadata = args->metadata;
return RES_OK;
}
+static res_T
+setup_lines(struct sln_tree* tree, const char* caller)
+{
+ size_t iline, nlines;
+ res_T res = RES_OK;
+ ASSERT(tree && caller);
+
+ SHTR(lines_view_get_size(tree->lines_view, &nlines));
+
+ res = darray_line_resize(&tree->lines, nlines);
+ if(res != RES_OK) {
+ log_err(tree->sln,
+ "%s: error allocating the list of precomputed line data.\n",
+ caller);
+ goto error;
+ }
+
+ FOR_EACH(iline, 0, nlines) {
+ line_setup(tree, iline);
+ }
+
+exit:
+ return res;
+error:
+ darray_line_clear(&tree->lines);
+ goto exit;
+}
+
static void
release_tree(ref_T* ref)
{
@@ -312,7 +344,9 @@ release_tree(ref_T* ref)
struct sln_device* sln = NULL;
ASSERT(ref);
sln = tree->sln;
+ if(tree->metadata) SHTR(isotope_metadata_ref_put(tree->metadata));
if(tree->lines_view) SHTR(lines_view_ref_put(tree->lines_view));
+ darray_line_release(&tree->lines);
darray_node_release(&tree->nodes);
darray_vertex_release(&tree->vertices);
MEM_RM(sln->allocator, tree);
@@ -345,17 +379,17 @@ sln_tree_create
ref_init(&tree->ref);
SLN(device_ref_get(sln));
tree->sln = sln;
+ darray_line_init(sln->allocator, &tree->lines);
darray_node_init(sln->allocator, &tree->nodes);
darray_vertex_init(sln->allocator, &tree->vertices);
- res = create_lines_view(sln, args, &tree->lines_view);
- if(res != RES_OK) goto error;
- res = setup_molecule_params(sln, args, tree->molecules_params);
- if(res != RES_OK) goto error;
- res = store_input_args(sln, args, tree);
- if(res != RES_OK) goto error;
- res = tree_build(tree, args);
- if(res != RES_OK) goto error;
+ #define CALL(Func) { if(RES_OK != (res = Func)) goto error; } (void)0
+ CALL(create_lines_view(sln, args, &tree->lines_view));
+ CALL(setup_molecule_params(sln, args, tree->molecules_params));
+ CALL(store_input_args(sln, args, tree));
+ CALL(setup_lines(tree, FUNC_NAME));
+ CALL(tree_build(tree, args));
+ #undef CALL
exit:
if(out_tree) *out_tree = tree;
diff --git a/src/sln_tree_c.h b/src/sln_tree_c.h
@@ -19,6 +19,9 @@
#ifndef SLN_TREE_C_H
#define SLN_TREE_C_H
+#include "sln.h"
+#include "sln_line.h"
+
#include <rsys/dynamic_array.h>
#include <rsys/ref_count.h>
@@ -42,6 +45,11 @@ struct vertex { /* 8 Bytes */
float ka;
};
+/* Generate the dynamic array of lines */
+#define DARRAY_DATA struct line
+#define DARRAY_NAME line
+#include <rsys/dynamic_array.h>
+
/* Generate the dynamic array of nodes */
#define DARRAY_DATA struct sln_node
#define DARRAY_NAME node
@@ -65,11 +73,18 @@ struct sln_tree {
* isotopes abundance) */
struct molecule_params molecules_params[SLN_MAX_MOLECULES_COUNT];
+ struct shtr_isotope_metadata* metadata;
struct shtr_lines_view* lines_view; /* Set of lines */
+
+ /* Store per line precomputed data */
+ struct darray_line lines;
+
struct darray_node nodes; /* Nodes used to partition the lines */
struct darray_vertex vertices; /* List of vertices */
- double wavenumber_range[2]; /* Spectral range */
+ double temperature; /* In Kelvin */
+ double pressure; /* In atm */
+ double wavenumber_range[2]; /* Spectral range [cm^-1] */
size_t max_nlines_per_leaf;
struct sln_device* sln;
@@ -81,4 +96,13 @@ tree_build
(struct sln_tree* tree,
const struct sln_tree_create_args* args);
+static INLINE const struct molecule_params*
+tree_get_molecule_params
+ (const struct sln_tree* tree,
+ const int32_t molecules_id)
+{
+ ASSERT(tree && molecules_id >= 0 && molecules_id < SLN_MAX_MOLECULES_COUNT);
+ return tree->molecules_params + molecules_id;
+}
+
#endif /* SLN_TREE_C_H */
