commit 22df719f66edfa9baba1fb2d8906b4b8d01c79bb
parent fb44c28d0176fa571d5c96bd754e4c27903c45a0
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 31 Jul 2018 16:07:06 +0200
Test the interpolation of the radiative properties
Diffstat:
8 files changed, 229 insertions(+), 5 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -100,11 +100,14 @@ if(NOT NO_TEST)
endif()
new_test(test_htgop)
+ build_test(test_htgop_fetch_radiative_properties)
build_test(test_htgop_load)
build_test(test_htgop_sample)
add_test(test_htgop_load test_htgop_load ${_etc_dst}/ecrad_opt_prop.txt)
add_test(test_htgop_sample test_htgop_sample ${_etc_dst}/ecrad_opt_prop.txt)
+ add_test(test_htgop_fetch_radiative_properties
+ test_htgop_fetch_radiative_properties ${_etc_dst}/ecrad_opt_prop.txt)
endif()
################################################################################
diff --git a/src/htgop.c b/src/htgop.c
@@ -275,11 +275,24 @@ load_stream(struct htgop* htgop, FILE* stream, const char* stream_name)
FOR_EACH(ilay, 0, nlays) {
double* x_h2o_tab = darray_double_data_get(&layers[ilay].x_h2o_tab);
CALL(cstr_to_double(read_line(&rdr), &x_h2o_tab[itab]));
-
+ if(x_h2o_tab[itab] < 0) {
+ log_err(htgop,
+ "%s:%lu: the tabulated water vapor molar fractions must be in [0, 1].\n",
+ rdr.name, rdr.iline);
+ res = RES_BAD_ARG;
+ goto error;
+ }
if(itab && x_h2o_tab[itab] < x_h2o_tab[itab-1]) {
log_err(htgop,
- "%s:%lu: The tabulated x H2O must be sorted in strict ascending "
- "order.\n", rdr.name, rdr.iline);
+ "%s:%lu: the tabulated water vapor molar fractions must be sorted in "
+ "strict ascending order.\n", rdr.name, rdr.iline);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ if(itab == tab_len-1 && x_h2o_tab[itab] != 1) {
+ log_err(htgop,
+ "%s:%lu: the tabulated water vapor molar fractions is not normalized.\n",
+ rdr.name, rdr.iline);
res = RES_BAD_ARG;
goto error;
}
diff --git a/src/htgop.h b/src/htgop.h
@@ -216,7 +216,7 @@ htgop_layer_sw_spectral_interval_get_tab
struct htgop_layer_sw_spectral_interval_tab* tab);
/*******************************************************************************
- * The following functions Retrieve the absorption/scattering coefficient in a
+ * The following functions retrieve the absorption/scattering coefficient in a
* layer, in short/long wave, at a given spectral interval for a specific
* quadrature point with respect to a submitted water vapor molar fraction
* x_h2o. x_h2o is compared to the values of tabulated water vapor molar
diff --git a/src/htgop_fetch_radiative_properties.h b/src/htgop_fetch_radiative_properties.h
@@ -83,7 +83,7 @@ CONCAT(CONCAT(CONCAT(htgop_layer_fetch_,LAYER_SPECINT_DOMAIN),_),LAYER_SPECINT_D
if(!k1 || !k2) { /* Linear interpolation */
k = k1 + (k2-k1)/(x2-x1) * (x_h2o-x1);
} else {
- const double alpha = (log(k2) - log(k1))/(log(x1) - log(x2));
+ const double alpha = (log(k2) - log(k1))/(log(x2) - log(x1));
const double beta = log(k1) - alpha * log(x1);
k = exp(alpha * log(x_h2o) + beta);
}
diff --git a/src/htgop_parse_layers_spectral_intervals_data.h b/src/htgop_parse_layers_spectral_intervals_data.h
@@ -85,6 +85,13 @@ CONCAT(parse_layers_spectral_intervals_, DATA)
(&layers[ilay].XDOMAIN(specints)) + ispecint;
XDATA(tab) = darray_dbllst_data_get(&layspecint->XDATA(tab)) + iquad;
CALL(cstr_to_double(read_line(rdr), &darray_double_data_get(XDATA(tab))[itab]));
+#if 0
+ if(darray_double_data_get(XDATA(tab))[itab] < 0) {
+ log_err(htgop, "The radiative properties cannot be negative.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+#endif
}
}
}
diff --git a/src/test_htgop_fetch_radiative_properties.c b/src/test_htgop_fetch_radiative_properties.c
@@ -0,0 +1,69 @@
+/* Copyright (C) |Meso|Star> 2018 (contact@meso-star.com)
+ *
+ * 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 "htgop.h"
+
+#include "htgop.h"
+#include "test_htgop_utils.h"
+
+#include <rsys/math.h>
+
+/* Generate the check_layer_fetch_lw_ka function */
+#define DOMAIN lw
+#define DATA ka
+#include "test_htgop_fetch_radiative_properties.h"
+/* Generate the check_layer_fetch_sw_ka function */
+#define DOMAIN sw
+#define DATA ka
+#include "test_htgop_fetch_radiative_properties.h"
+/* Generate the check_layer_fetch_sw_ks function */
+#define DOMAIN sw
+#define DATA ks
+#include "test_htgop_fetch_radiative_properties.h"
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct htgop* htgop;
+ size_t i, n;
+
+ if(argc < 2) {
+ fprintf(stderr, "Usage: %s FILENAME\n", argv[0]);
+ return 1;
+ }
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(htgop_create(NULL, &allocator, 1, &htgop) == RES_OK);
+ CHK(htgop_load(htgop, argv[1]) == RES_OK);
+
+ CHK(htgop_get_layers_count(htgop, &n) == RES_OK);
+ CHK(n > 0);
+
+ FOR_EACH(i, 0, n) {
+ struct htgop_layer layer;
+ CHK(htgop_get_layer(htgop, i, &layer) == RES_OK);
+ check_layer_fetch_lw_ka(htgop, &layer);
+ check_layer_fetch_sw_ka(htgop, &layer);
+ check_layer_fetch_sw_ks(htgop, &layer);
+ }
+
+ CHK(htgop_ref_put(htgop) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_htgop_fetch_radiative_properties.h b/src/test_htgop_fetch_radiative_properties.h
@@ -0,0 +1,128 @@
+/* Copyright (C) |Meso|Star> 2018 (contact@meso-star.com)
+ *
+ * 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 "htgop.h"
+#include <math.h>
+
+#if !defined(DOMAIN) || !defined(DATA)
+ #error "Missing the <DATA|DOMAIN> macro."
+#endif
+
+/* Helper macros */
+#define GET_NSPECINTS \
+ CONCAT(CONCAT(htgop_get_, DOMAIN), _spectral_intervals_count)
+#define GET_SPECINT \
+ CONCAT(CONCAT(htgop_get_, DOMAIN), _spectral_interval)
+#define LAY_SPECINT \
+ CONCAT(CONCAT(htgop_layer_, DOMAIN),_spectral_interval)
+#define LAY_SPECINT_TAB \
+ CONCAT(CONCAT(htgop_layer_, DOMAIN),_spectral_interval_tab)
+#define GET_LAY_SPECINT \
+ CONCAT(CONCAT(htgop_layer_get_, DOMAIN), _spectral_interval)
+#define GET_LAY_SPECINT_TAB \
+ CONCAT(CONCAT(htgop_layer_, DOMAIN), _spectral_interval_get_tab)
+#define FETCH_K \
+ CONCAT(CONCAT(CONCAT(htgop_layer_fetch_, DOMAIN),_), DATA)
+#define K_TAB CONCAT(DATA, _tab)
+
+static void
+CONCAT(CONCAT(CONCAT(check_layer_fetch_, DOMAIN),_), DATA)
+ (const struct htgop* htgop,
+ const struct htgop_layer* layer)
+{
+ struct htgop_spectral_interval specint;
+ size_t ispecint;
+ size_t nspecints;
+ size_t quad_len;
+ double x_h2o;
+ double k;
+
+ CHK(htgop && layer);
+ CHK(GET_NSPECINTS(htgop, &nspecints) == RES_OK);
+ CHK(nspecints);
+
+ CHK(GET_SPECINT(htgop, 0, &specint) == RES_OK);
+ CHK(specint.quadrature_length);
+ CHK(layer->tab_length);
+
+ quad_len = specint.quadrature_length;
+
+ x_h2o = layer->x_h2o_tab[0];
+ CHK(FETCH_K(NULL, 0, 0, x_h2o, &k) == RES_BAD_ARG);
+ CHK(FETCH_K(layer, nspecints, 0, x_h2o, &k) == RES_BAD_ARG);
+ CHK(FETCH_K(layer, 0, quad_len, x_h2o, &k) == RES_BAD_ARG);
+ CHK(FETCH_K(layer, 0, 0, 1.0000001, &k) == RES_BAD_ARG);
+ CHK(FETCH_K(layer, 0, 0, x_h2o, NULL) == RES_BAD_ARG);
+
+ FOR_EACH(ispecint, 0, nspecints) {
+ struct LAY_SPECINT lay_specint;
+ struct LAY_SPECINT_TAB lay_tab;
+ size_t iquad;
+
+ CHK(GET_SPECINT(htgop, ispecint, &specint) == RES_OK);
+ CHK(specint.quadrature_length);
+
+ FOR_EACH(iquad, 0, specint.quadrature_length) {
+ size_t itab;
+
+ CHK(GET_LAY_SPECINT(layer, ispecint, &lay_specint) == RES_OK);
+ CHK(GET_LAY_SPECINT_TAB(&lay_specint, iquad, &lay_tab) == RES_OK);
+ CHK(lay_tab.tab_length == layer->tab_length);
+
+ /* Check boundaries */
+ FOR_EACH(itab, 0, layer->tab_length) {
+ const size_t N = 10;
+ size_t i;
+
+ x_h2o = layer->x_h2o_tab[itab];
+ CHK(FETCH_K(layer, ispecint, iquad, x_h2o, &k) == RES_OK);
+ CHK(eq_eps(k, lay_tab.K_TAB[itab], 1.e-6));
+
+ FOR_EACH(i, 0, N) {
+ const double r = rand_canonic();
+ const double x1 = itab ? layer->x_h2o_tab[itab-1] : 0;
+ const double k1 = itab ? lay_tab.K_TAB[itab-1] : 0;
+ const double x2 = layer->x_h2o_tab[itab];
+ const double k2 = lay_tab.K_TAB[itab];
+ double ref;
+
+ x_h2o = x1 + (x2 - x1) * r;
+ if(!itab || !k1 || !k2) {
+ ref = k1 + (k2-k1) / (x2 - x1) * x_h2o;
+ } else {
+ const double alpha = (log(k2) - log(k1)) / (log(x2) - log(x1));
+ const double beta = log(k1) - alpha*log(x1);
+ ref = exp(alpha*log(x_h2o) + beta);
+ }
+
+ CHK(FETCH_K(layer, ispecint, iquad, x_h2o, &k) == RES_OK);
+ CHK(eq_eps(k, ref, 1.e-6));
+ }
+ }
+ }
+ }
+}
+
+#undef GET_NSPECINTS
+#undef GET_SPECINT
+#undef LAY_SPECINT
+#undef LAY_SPECINT_TAB
+#undef GET_LAY_SPECINT
+#undef GET_LAY_SPECINT_TAB
+#undef FETCH_K
+#undef K_TAB
+#undef DOMAIN
+#undef DATA
+
diff --git a/src/test_htgop_load.c b/src/test_htgop_load.c
@@ -288,6 +288,8 @@ main(int argc, char** argv)
FOR_EACH(iquad, 0, specint.quadrature_length) {
FOR_EACH(ilay, 0, nlays) {
CHK(htgop_get_layer(htgop, ilay, &lay) == RES_OK);
+ CHK(lay.lw_spectral_intervals_count == nspecints_lw);
+ CHK(lay.sw_spectral_intervals_count == nspecints_sw);
CHK(htgop_layer_get_lw_spectral_interval(&lay, ispecint, &lw_specint) == RES_OK);
CHK(lw_specint.quadrature_length == specint.quadrature_length);
CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK);
@@ -296,6 +298,8 @@ main(int argc, char** argv)
FOR_EACH(itab, 0, tab_len) {
FOR_EACH(ilay, 0, nlays) {
CHK(htgop_get_layer(htgop, ilay, &lay) == RES_OK);
+ CHK(lay.lw_spectral_intervals_count == nspecints_lw);
+ CHK(lay.sw_spectral_intervals_count == nspecints_sw);
CHK(htgop_layer_get_lw_spectral_interval(&lay, ispecint, &lw_specint) == RES_OK);
CHK(htgop_layer_lw_spectral_interval_get_tab(&lw_specint, iquad, &lw_tab) == RES_OK);
CHK(lw_tab.tab_length == tab_len);
