commit 1bb634ce9ac26608fd1a6a4fcd4d448efd5130db
parent abdbd1f6aa716905a4eba1d26dadb8ea2e5707b7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 25 May 2022 17:46:49 +0200
Generate a mesh that is an upper limit of the line
Diffstat:
3 files changed, 184 insertions(+), 29 deletions(-)
diff --git a/src/sln_line.c b/src/sln_line.c
@@ -16,6 +16,8 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
+#define _POSIX_C_SOURCE 200112L /* nextafterf support */
+
#include "sln_device_c.h"
#include "sln_line.h"
#include "sln_log.h"
@@ -25,10 +27,13 @@
#include <lblu.h>
#include <star/shtr.h>
+#include <rsys/algorithm.h>
#include <rsys/cstr.h>
#include <rsys/dynamic_array_double.h>
#include <rsys/math.h>
+#include <math.h> /* nextafterf */
+
#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 */
@@ -240,6 +245,8 @@ eval_mesh
const enum sln_line_profile line_profile,
struct darray_double* values)
{
+ const double* nu = NULL;
+ double* ha = NULL;
const struct shtr_line* shtr_line = NULL;
size_t ivertex, nvertices;
double range[2];
@@ -267,12 +274,22 @@ eval_mesh
range[0] = MMIN(shtr_line->wavenumber, mixture->wavenumber_range[0]);
range[1] = shtr_line->wavenumber + cutoff;
+ nu = darray_double_cdata_get(wavenumbers);
+ ha = darray_double_data_get(values);
FOR_EACH(ivertex, 0, nvertices) {
- const double nu = darray_double_cdata_get(wavenumbers)[ivertex];
- if(nu < range[0]) continue;
- if(nu > range[1]) break;
- darray_double_data_get(values)[ivertex] =
- line_value(mixture, iline, line_profile, nu);
+ 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(mixture, iline, line_profile, nu_curr);
+ }
}
exit:
@@ -283,6 +300,79 @@ error:
goto exit;
}
+static void
+snap_mesh_to_upper_bound
+ (const struct darray_double* wavenumbers,
+ struct darray_double* values)
+{
+ double* ha = NULL;
+ size_t ivertex, nvertices;
+ size_t ivertex_1st;
+ ASSERT(wavenumbers && values);
+ ASSERT(darray_double_size_get(wavenumbers) == darray_double_size_get(values));
+ (void)wavenumbers;
+
+ ha = darray_double_data_get(values);
+ nvertices = darray_double_size_get(wavenumbers);
+
+ /* Search for the first (non clipped) vertex */
+ FOR_EACH(ivertex_1st, 0, nvertices) {
+ if(ha[ivertex_1st] > 0) break;
+ }
+ ASSERT(ivertex_1st < nvertices);
+
+ /* Ensure that the stored vertex value is an exclusive upper bound of the
+ * original value. We do this by storing a value in single precision that is
+ * strictly greater than its encoding in double precision */
+ if(ha[ivertex_1st] != (float)ha[ivertex_1st]) {
+ ha[ivertex_1st] = nextafterf((float)ha[ivertex_1st], FLT_MAX);
+ }
+
+ /* We have meshed the upper half of the line which is a strictly decreasing
+ * function. To ensure that the mesh is an upper limit of this function,
+ * simply align the value of each vertex with the value of the preceding
+ * vertex */
+ FOR_EACH_REVERSE(ivertex, nvertices-1, ivertex_1st) {
+ if(ha[ivertex] < 0) continue; /* The vertex is clipped */
+ ha[ivertex] = ha[ivertex-1];
+ }
+}
+
+static INLINE int
+cmp_dbl(const void* a, const void* b)
+{
+ const double key = *((const double*)a);
+ const double item = *((const double*)b);
+ if(key < item) return -1;
+ if(key > item) return +1;
+ return 0;
+}
+
+/* Return the value of the vertex whose wavenumber is greater than 'nu' */
+static INLINE double
+next_vertex_value
+ (const double nu,
+ const struct darray_double* wavenumbers,
+ const struct darray_double* values)
+{
+ const double* wnum = NULL;
+ size_t ivertex = 0;
+ ASSERT(wavenumbers && values);
+
+ wnum = search_lower_bound
+ (&nu,
+ darray_double_cdata_get(wavenumbers),
+ darray_double_size_get(wavenumbers),
+ sizeof(double),
+ cmp_dbl);
+ ASSERT(wnum); /* It necessary exists */
+
+ ivertex = (size_t)(wnum - darray_double_cdata_get(wavenumbers));
+ ASSERT(ivertex < darray_double_size_get(values));
+
+ return darray_double_cdata_get(values)[ivertex];
+}
+
/* Append the line mesh into the vertices array */
static res_T
save_line_mesh
@@ -292,17 +382,22 @@ save_line_mesh
const double spectral_range[2],
const struct darray_double* wavenumbers,
const struct darray_double* values,
+ const enum mesh_type mesh_type,
size_t vertices_range[2]) /* Range into which the line vertices are saved */
{
const struct line* line = NULL;
+ const double* wnums = NULL;
+ const double* vals = 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_range);
ASSERT(darray_double_size_get(wavenumbers) == darray_double_size_get(values));
+ ASSERT((unsigned)mesh_type < MESH_TYPES_COUNT__);
ASSERT(spectral_range && spectral_range[0] <= spectral_range[1]);
line = darray_line_cdata_get(&tree->mixture->lines) + iline;
@@ -322,62 +417,104 @@ save_line_mesh
goto error;
}
+ wnums = darray_double_cdata_get(wavenumbers);
+ vals = darray_double_cdata_get(values);
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;
- vtx->wavenumber = (float)spectral_range[0];
- vtx->ka = (float)line_value
- (tree->mixture, iline, tree->line_profile, vtx->wavenumber);
+ struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices)+i;
+ const double nu = spectral_range[0];
+ double ha = 0;
+
+ switch(mesh_type) {
+ case MESH_UPPER_BOUND:
+ ha = nu > line->wavenumber
+ ? next_vertex_value(nu, wavenumbers, values)
+ : next_vertex_value(MIRROR(nu), wavenumbers, values);
+ break;
+
+ case MESH_USUAL:
+ ha = line_value(tree->mixture, iline, tree->line_profile, 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;
- const double ka = darray_double_cdata_get(values)[ivertex];
- const double nu = /* Mirror the wavenumber */
- 2*line->wavenumber - darray_double_cdata_get(wavenumbers)[ivertex];
+ struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices)+i;
+ const double nu = 2*line->wavenumber - wnums[ivertex]; /* Mirror */
+ 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;
+ 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)ka;
+ 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;
- const double ka = darray_double_cdata_get(values)[ivertex];
- const double nu = darray_double_cdata_get(wavenumbers)[ivertex];
+ 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;
+ if(nu >= spectral_range[1]) break;
vtx->wavenumber = (float)nu;
- vtx->ka = (float)ka;
+ 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;
- vtx->wavenumber = (float)spectral_range[1];
- vtx->ka = (float)line_value
- (tree->mixture, iline, tree->line_profile, vtx->wavenumber);
+ struct sln_vertex* vtx = darray_vertex_data_get(&tree->vertices)+i;
+ const double nu = spectral_range[1];
+ double ha = 0;
+
+ switch(mesh_type) {
+ case MESH_UPPER_BOUND:
+ ha = nu > line->wavenumber
+ ? next_vertex_value(nu, wavenumbers, values)
+ : next_vertex_value(MIRROR(nu), wavenumbers, values);
+ break;
+
+ case MESH_USUAL:
+ ha = line_value(tree->mixture, iline, tree->line_profile, 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. */
@@ -476,6 +613,7 @@ line_mesh
(struct sln_tree* tree,
const size_t iline,
const size_t nvertices_hint,
+ const enum mesh_type mesh_type,
size_t vertices_range[2])
{
struct darray_double values; /* List of evaluated values */
@@ -487,6 +625,7 @@ line_mesh
res_T res = RES_OK;
ASSERT(tree && nvertices_hint);
ASSERT(iline < darray_line_size_get(&tree->mixture->lines));
+ ASSERT((unsigned)mesh_type < MESH_TYPES_COUNT__);
line = darray_line_cdata_get(&tree->mixture->lines) + iline;
SHTR(line_view_get_line(tree->mixture->line_view, iline, &shtr_line));
@@ -513,8 +652,16 @@ line_mesh
* wavenumbers */
eval_mesh(tree->mixture, iline, &wavenumbers, tree->line_profile, &values);
+ switch(mesh_type) {
+ case MESH_UPPER_BOUND:
+ snap_mesh_to_upper_bound(&wavenumbers, &values);
+ break;
+ case MESH_USUAL: /* Do nothing */ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+
res = save_line_mesh(tree, iline, mol_params, tree->mixture->wavenumber_range,
- &wavenumbers, &values, vertices_range);
+ &wavenumbers, &values, mesh_type, vertices_range);
if(res != RES_OK) goto error;
exit:
diff --git a/src/sln_line.h b/src/sln_line.h
@@ -32,6 +32,12 @@ struct line {
double gamma_l; /* Lorentz half width */
};
+enum mesh_type {
+ MESH_UPPER_BOUND,
+ MESH_USUAL,
+ MESH_TYPES_COUNT__
+};
+
/* Forward declaration */
struct sln_mixture;
struct sln_tree;
@@ -63,7 +69,8 @@ line_mesh
(struct sln_tree* tree,
const size_t iline,
const size_t nvertices_hint,
- /* Range of generated line vertices. Upper bound is exclusive */
+ const enum mesh_type mesh_type,
+ /* Range of generated line vertices. Upper bound is inclusive */
size_t vertices_range[2]);
#endif /* SLN_LINE_H */
diff --git a/src/sln_tree_build.c b/src/sln_tree_build.c
@@ -55,7 +55,8 @@ build_leaf_polyline
ASSERT(leaf->range[0] == leaf->range[1]);
/* Line meshing */
- res = line_mesh(tree, leaf->range[0], args->nvertices_hint, vertices_range);
+ res = line_mesh(tree, leaf->range[0], args->nvertices_hint, MESH_USUAL,
+ vertices_range);
if(res != RES_OK) goto error;
/* Decimate the line mesh */
