commit 551b5ce63b058bdb4e2acebe09bd1fa94bbda6d6
parent 860df195e90c9592215f792d51e1d75d345fb755
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 27 Mar 2024 15:07:07 +0100
Fix use of unknown temperature
Since commit 75ae9b8, the unknown temperature is defined by a constant
which is currently set to "Not a Number". We use NaN to finely check
that no unknown temperature is used anywhere, since a NaN in the
calculation will lead to a NaN as the final result. This is a very
strict constraint, as it forbids inhibiting an unknown temperature by
multiplying it by 0. And this commit solves a number of cases of this
kind which, on the face of it, are not problematic. Nevertheless, even
if the use of NaN seems very strict, the effort to avoid using any
unknown temperature has highlighted several bugs or lack of verification
that this commit corrects.
In addition, some constants were not correctly set to unknown
temperatures, such as the temperature range of default scene creation
arguments.
During this review and debug, we also corrected a concurrent write to a
shared variable when calculating flux at a probe boundary. Strictly
speaking, this correction is not related to the use of an unknown
temperature, but it occurred in the middle of a rewrite of the same
piece of code, for which it seems artificial and excessive to rework
separately to make a clean commit.
Diffstat:
10 files changed, 180 insertions(+), 53 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -53,7 +53,7 @@
/* Syntactic sugar used to define whether a temperature is known or not */
#define SDIS_TEMPERATURE_NONE NaN /* Unknown temperature */
#define SDIS_TEMPERATURE_IS_KNOWN(Temp) (!IS_NaN(Temp))
-#define SDIS_TEMPERATURE_IS_UNKNOWN(Temp) (!SDIS_TEMPERATURE_IS_KNOWN(Temp))
+#define SDIS_TEMPERATURE_IS_UNKNOWN(Temp) (IS_NaN(Temp))
/* Forward declaration of external opaque data types */
struct logger;
@@ -454,7 +454,10 @@ struct sdis_ambient_radiative_temperature {
double temperature; /* In Kelvin */
double reference; /* Used to linearise the radiative transfer */
};
-#define SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__ {-1, -1}
+#define SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__ { \
+ SDIS_TEMPERATURE_NONE, \
+ SDIS_TEMPERATURE_NONE \
+}
static const struct sdis_ambient_radiative_temperature
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL =
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__;
@@ -490,7 +493,7 @@ struct sdis_scene_create_args {
0, /* #vertices */ \
1.0, /* #Floating point to meter scale factor */ \
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__,/* Ambient radiative temperature */\
- {0.0, -1.0}, /* Temperature range */ \
+ {SDIS_TEMPERATURE_NONE, SDIS_TEMPERATURE_NONE}, /* Temperature range */ \
NULL /* source */ \
}
static const struct sdis_scene_create_args SDIS_SCENE_CREATE_ARGS_DEFAULT =
diff --git a/src/sdis_heat_path_boundary_Xd_c.h b/src/sdis_heat_path_boundary_Xd_c.h
@@ -992,4 +992,40 @@ error:
goto exit;
}
+res_T
+XD(check_Tref)
+ (const struct sdis_scene* scn,
+ const double pos[DIM],
+ const double Tref,
+ const char* func_name)
+{
+ ASSERT(scn && pos && func_name);
+
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(Tref)) {
+ log_err(scn->dev,
+ "%s: invalid reference temperature at the position `"FORMAT_VECX"'.\n",
+ func_name, SPLITX(pos));
+ return RES_BAD_OP_IRRECOVERABLE;
+ }
+
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(scn->tmin)
+ || SDIS_TEMPERATURE_IS_UNKNOWN(scn->tmax)) {
+ log_err(scn->dev,
+ "%s: invalid scene temperature range. "
+ "At least one boundary is unknown.\n",
+ func_name);
+ return RES_BAD_OP_IRRECOVERABLE;
+ }
+
+ if(Tref > scn->tmax) {
+ log_err(scn->dev,
+ "%s: invalid maximum temperature `%gK'. The reference temperature `%gK' "
+ "at the position `"FORMAT_VECX"' is greater than this temperature.\n",
+ func_name, scn->tmax, Tref, SPLITX(pos));
+ return RES_BAD_OP_IRRECOVERABLE;
+ }
+
+ return RES_OK;
+}
+
#include "sdis_Xd_end.h"
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h
@@ -28,32 +28,6 @@
* Helper functions
******************************************************************************/
static INLINE res_T
-XD(check_Tref)
- (const struct sdis_scene* scn,
- const double pos[DIM],
- const double Tref,
- const char* func_name)
-{
- ASSERT(scn && pos && func_name);
-
- if(SDIS_TEMPERATURE_IS_UNKNOWN(Tref)) {
- log_err(scn->dev,
- "%s: invalid reference temperature at the position `"FORMAT_VECX"'.\n",
- func_name, SPLITX(pos));
- return RES_BAD_OP_IRRECOVERABLE;
- }
- if(Tref > scn->tmax) {
- log_err(scn->dev,
- "%s: invalid maximum temperature `%gK'. The reference temperature `%gK' "
- "at the position `"FORMAT_VECX"' is greater than this temperature.\n",
- func_name, scn->tmax, Tref, SPLITX(pos));
- return RES_BAD_OP_IRRECOVERABLE;
- }
-
- return RES_OK;
-}
-
-static INLINE res_T
XD(rwalk_get_Tref)
(const struct sdis_scene* scn,
const struct XD(rwalk)* rwalk,
@@ -213,7 +187,13 @@ XD(solid_fluid_boundary_picard1_path)
/* Compute the convective, conductive and the upper bound radiative coef */
h_conv = interface_get_convection_coef(interf, frag);
h_cond = lambda / delta_m;
- h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * ctx->That3 * epsilon;
+ if(epsilon <= 0) {
+ h_radi_hat = 0; /* No radiative transfert */
+ } else {
+ res = scene_check_temperature_range(scn);
+ if(res != RES_OK) { res = RES_BAD_OP_IRRECOVERABLE; goto error; }
+ h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * ctx->That3 * epsilon;
+ }
/* Compute a global upper bound coefficient */
h_hat = h_conv + h_cond + h_radi_hat;
@@ -304,6 +284,11 @@ XD(solid_fluid_boundary_picard1_path)
/* Get the Tref at the end of the candidate radiative path */
res = XD(rwalk_get_Tref)(scn, &rwalk_s, &T_s, &Tref_s);
if(res != RES_OK) goto error;
+
+ /* The reference temperatures must be known, as this is a radiative path.
+ * If this is not the case, an error should be reported before this point.
+ * Hence these assertions to detect unexpected behavior */
+ ASSERT(SDIS_TEMPERATURE_IS_KNOWN(Tref));
ASSERT(SDIS_TEMPERATURE_IS_KNOWN(Tref_s));
h_radi = BOLTZMANN_CONSTANT * epsilon *
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h
@@ -232,7 +232,15 @@ XD(solid_fluid_boundary_picardN_path)
/* Compute the convective, conductive and the upper bound radiative coef */
h_conv = interface_get_convection_coef(interf, frag);
h_cond = lambda / (delta * scn->fp_to_meter);
- h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * That3 * epsilon;
+ h_radi_hat = epsilon > 0 ? 4.0 * BOLTZMANN_CONSTANT * That3 * epsilon : 0;
+
+ if(epsilon <= 0) {
+ h_radi_hat = 0; /* No radiative transfert */
+ } else {
+ res = scene_check_temperature_range(scn);
+ if(res != RES_OK) { res = RES_BAD_OP_IRRECOVERABLE; goto error; }
+ h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * That3 * epsilon;
+ }
/* Compute a global upper bound coefficient */
h_hat = h_conv + h_cond + h_radi_hat;
diff --git a/src/sdis_heat_path_boundary_c.h b/src/sdis_heat_path_boundary_c.h
@@ -221,6 +221,23 @@ handle_external_net_flux_3d
struct temperature_3d* T);
/*******************************************************************************
+ * Miscellaneous functions
+ ******************************************************************************/
+extern LOCAL_SYM res_T
+check_Tref_2d
+ (const struct sdis_scene* scn,
+ const double pos[2],
+ const double Tref,
+ const char* call_func_name);
+
+extern LOCAL_SYM res_T
+check_Tref_3d
+ (const struct sdis_scene* scn,
+ const double pos[3],
+ const double Tref,
+ const char* call_func_name);
+
+/*******************************************************************************
* Boundary sub-paths
******************************************************************************/
extern LOCAL_SYM res_T
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -578,3 +578,41 @@ exit:
error:
goto exit;
}
+
+res_T
+scene_check_temperature_range(const struct sdis_scene* scn)
+{
+ res_T res = RES_OK;
+ ASSERT(scn);
+
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(scn->tmin)) {
+ log_err(scn->dev,
+ "%s the defined minimum temperature is unknown "
+ "when it is expected to be known.\n",
+ FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(scn->tmax)) {
+ log_err(scn->dev,
+ "%s the defined maximum temperature is unknown "
+ "when it is expected to be known.\n",
+ FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(scn->tmin > scn->tmax) {
+ log_err(scn->dev,
+ "%s: defined temperature range degenerated -- [%g, %g] K\n",
+ FUNC_NAME, scn->tmin, scn->tmax);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -290,6 +290,14 @@ extern LOCAL_SYM res_T
scene_check_dimensionality_3d
(const struct sdis_scene* scn);
+/* Check that the temperature range of the scene is well defined, i.e. that the
+ * minimum and maximum temperatures are known and that they define a valid
+ * range. If this is not the case, the function displays an error message and
+ * returns RES_BAD_ARG */
+extern LOCAL_SYM res_T
+scene_check_temperature_range
+ (const struct sdis_scene* scn);
+
static INLINE void
scene_get_enclosure_ids
(const struct sdis_scene* scn,
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -23,6 +23,7 @@
#include "sdis_misc.h"
#include "sdis_realisation.h"
#include "sdis_scene_c.h"
+#include "sdis_heat_path_boundary_c.h" /* check_Tref_<2d|3d> */
#include <rsys/clock_time.h>
#include <star/ssp.h>
@@ -784,9 +785,16 @@ XD(solve_boundary_flux)
/* Fetch interface parameters */
epsilon = interface_side_get_emissivity(interf, &frag);
- Tref = interface_side_get_reference_temperature(interf, &frag);
hc = interface_get_convection_coef(interf, &frag);
- hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
+ Tref = interface_side_get_reference_temperature(interf, &frag);
+ if(epsilon <= 0) {
+ hr = 0;
+ } else {
+ res_local = XD(check_Tref)(scn, frag.P, Tref, FUNC_NAME);
+ if(res_local != RES_OK) { ATOMIC_SET(&res, &res_local); continue; }
+ hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
+ }
+
frag.side = solid_side;
imposed_flux = interface_side_get_flux(interf, &frag);
imposed_temp = interface_side_get_temperature(interf, &frag);
@@ -827,10 +835,9 @@ XD(solve_boundary_flux)
} else if(res_simul == RES_OK) { /* Update accumulators */
const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
/* Convective flux from fluid to solid */
- const double w_conv = hc * (result.Tfluid - result.Tboundary);
+ const double w_conv = hc > 0 ? hc * (result.Tfluid - result.Tboundary) : 0;
/* Radiative flux from ambient to solid */
- const double w_rad = SDIS_TEMPERATURE_IS_UNKNOWN(result.Tradiative) ?
- 0 : hr * (result.Tradiative - result.Tboundary);
+ const double w_rad = hr > 0 ? hr * (result.Tradiative - result.Tboundary) : 0;
/* Imposed flux that goes _into_ the solid */
const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0;
const double w_total = w_conv + w_rad + w_imp;
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -23,6 +23,7 @@
#include "sdis_misc.h"
#include "sdis_realisation.h"
#include "sdis_scene_c.h"
+#include "sdis_heat_path_boundary_c.h" /* check_Tref_<2d|3d> */
#include <rsys/clock_time.h>
#include <star/ssp.h>
@@ -852,6 +853,7 @@ XD(solve_probe_boundary_flux)
BOUNDARY_FLUX_REALISATION_ARGS_NULL;
struct time t0, t1;
const int ithread = omp_get_thread_num();
+ struct sdis_interface_fragment frag_local = frag;
struct ssp_rng* rng = per_thread_rng[ithread];
struct accum* acc_temp = &per_thread_acc_tp[ithread];
struct accum* acc_time = &per_thread_acc_ti[ithread];
@@ -862,9 +864,10 @@ XD(solve_probe_boundary_flux)
double time, epsilon, hc, hr, imposed_flux, imposed_temp;
int flux_mask = 0;
struct bound_flux_result result = BOUND_FLUX_RESULT_NULL__;
- double Tref = -1;
+ double Tref = SDIS_TEMPERATURE_NONE;
size_t n;
int pcent;
+ res_T res_local = RES_OK;
res_T res_simul = RES_OK;
if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
@@ -875,15 +878,22 @@ XD(solve_probe_boundary_flux)
time = sample_time(rng, args->time_range);
/* Compute hr and hc */
- frag.time = time;
- frag.side = fluid_side;
- epsilon = interface_side_get_emissivity(interf, &frag);
- Tref = interface_side_get_reference_temperature(interf, &frag);
- hc = interface_get_convection_coef(interf, &frag);
- hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
- frag.side = solid_side;
- imposed_flux = interface_side_get_flux(interf, &frag);
- imposed_temp = interface_side_get_temperature(interf, &frag);
+ frag_local.time = time;
+ frag_local.side = fluid_side;
+ epsilon = interface_side_get_emissivity(interf, &frag_local);
+ Tref = interface_side_get_reference_temperature(interf, &frag_local);
+ hc = interface_get_convection_coef(interf, &frag_local);
+ if(epsilon <= 0) {
+ hr = 0;
+ } else {
+ res_local = XD(check_Tref)(scn, frag_local.P, Tref, FUNC_NAME);
+ if(res_local != RES_OK) { ATOMIC_SET(&res, &res_local); continue; }
+ hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon;
+ }
+
+ frag_local.side = solid_side;
+ imposed_flux = interface_side_get_flux(interf, &frag_local);
+ imposed_temp = interface_side_get_temperature(interf, &frag_local);
if(SDIS_TEMPERATURE_IS_KNOWN(imposed_temp)) {
/* Flux computation on T boundaries is not supported yet */
log_err(scn->dev,
@@ -922,10 +932,9 @@ XD(solve_probe_boundary_flux)
} else if(res_simul == RES_OK) { /* Update accumulators */
const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
/* Convective flux from fluid to solid */
- const double w_conv = hc * (result.Tfluid - result.Tboundary);
+ const double w_conv = hc > 0 ? hc * (result.Tfluid - result.Tboundary) : 0;
/* Radiative flux from ambient to solid */
- const double w_rad = SDIS_TEMPERATURE_IS_UNKNOWN(result.Tradiative) ?
- 0 : hr * (result.Tradiative - result.Tboundary);
+ const double w_rad = hr > 0 ? hr * (result.Tradiative - result.Tboundary) : 0;
/* Imposed flux that goes _into_ the solid */
const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0;
/* Total flux */
diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c
@@ -383,8 +383,16 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
BA(sdis_scene_get_ambient_radiative_temperature(scn, NULL));
BA(sdis_scene_get_ambient_radiative_temperature(NULL, &trad));
OK(sdis_scene_get_ambient_radiative_temperature(scn, &trad));
- CHK(trad.temperature == SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.temperature);
- CHK(trad.reference == SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.reference);
+ if(SDIS_TEMPERATURE_IS_KNOWN(trad.temperature)) {
+ CHK(trad.temperature == SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.temperature);
+ } else {
+ CHK(SDIS_TEMPERATURE_IS_UNKNOWN(SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.temperature));
+ }
+ if(SDIS_TEMPERATURE_IS_KNOWN(trad.reference)) {
+ CHK(trad.reference == SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.reference);
+ } else {
+ CHK(SDIS_TEMPERATURE_IS_UNKNOWN(SDIS_SCENE_CREATE_ARGS_DEFAULT.trad.reference));
+ }
trad.temperature = 100;
trad.reference = 110;
@@ -399,8 +407,16 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
BA(sdis_scene_get_temperature_range(scn, NULL));
BA(sdis_scene_get_temperature_range(NULL, t_range));
OK(sdis_scene_get_temperature_range(scn, t_range));
- CHK(t_range[0] == SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[0]);
- CHK(t_range[1] == SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[1]);
+ if(SDIS_TEMPERATURE_IS_KNOWN(t_range[0])) {
+ CHK(t_range[0] == SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[0]);
+ } else {
+ CHK(SDIS_TEMPERATURE_IS_UNKNOWN(SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[0]));
+ }
+ if(SDIS_TEMPERATURE_IS_KNOWN(t_range[1])) {
+ CHK(t_range[1] == SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[1]);
+ } else {
+ CHK(SDIS_TEMPERATURE_IS_UNKNOWN(SDIS_SCENE_CREATE_ARGS_DEFAULT.t_range[0]));
+ }
t_range[0] = 1;
t_range[1] = 100;
