commit 22f2c11ea95188399a306999b699fc448577115f
parent 51284033a44ed852765399fc26eda1efce1bd379
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 11 Sep 2024 11:56:40 +0200
Merge branch 'release_0.16.1'
Diffstat:
12 files changed, 244 insertions(+), 121 deletions(-)
diff --git a/README.md b/README.md
@@ -158,6 +158,17 @@ Edit config.mk as needed, then run:
## Release notes
+
+### Version 0.16.1
+
+- Corrected net flux calculation on surfaces with several Robin boundary
+ conditions. Depending on the configuration, such a calculation could
+ be impossible and return an error.
+- Improve the performance of external flux calculations by avoiding the
+ need to calculate the contribution of external flux for surfaces with
+ zero emissivity.
+- Mitigate numerical errors when sampling radiative paths.
+
### Version 0.16
#### Add support for custom sampling of solid paths
diff --git a/config.mk b/config.mk
@@ -1,6 +1,6 @@
VERSION_MAJOR = 0
VERSION_MINOR = 16
-VERSION_PATCH = 0
+VERSION_PATCH = 1
VERSION = $(VERSION_MAJOR).$(VERSION_MINOR).$(VERSION_PATCH)
PREFIX = /usr/local
diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h
@@ -25,58 +25,6 @@
#include "sdis_Xd_begin.h"
/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-/* This function checks whether the random walk is on a boundary and, if so,
- * verifies that the temperature of the medium attached to the interface is
- * known. This medium can be different from the medium of the enclosure. Indeed,
- * the enclosure can contain several media used to set the temperatures of
- * several boundary conditions. Hence this function, which queries the medium on
- * the trajectory coming from a boundary */
-static res_T
-XD(handle_known_medium_temperature)
- (struct sdis_scene* scn,
- struct rwalk_context* ctx,
- struct rwalk* rwalk,
- struct temperature* T)
-{
- struct sdis_interface* interf = NULL;
- struct sdis_medium* mdm = NULL;
- double temperature = SDIS_TEMPERATURE_NONE;
- res_T res = RES_OK;
- ASSERT(scn && ctx && rwalk && T);
-
- /* Not at an interface */
- if(SXD_HIT_NONE(&rwalk->XD(hit))) return RES_OK; /* Nothing to do */
-
- interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id);
- mdm = rwalk->hit_side==SDIS_FRONT ? interf->medium_front: interf->medium_back;
-
- temperature = medium_get_temperature(mdm, &rwalk->vtx);
-
- /* Check if the temperature is known */
- if(SDIS_TEMPERATURE_IS_UNKNOWN(temperature)) goto exit;
-
- T->value += temperature;
- T->done = 1;
-
- if(ctx->green_path) {
- res = green_path_set_limit_vertex
- (ctx->green_path, mdm, &rwalk->vtx, rwalk->elapsed_time);
- if(res != RES_OK) goto error;
- }
-
- if(ctx->heat_path) {
- heat_path_get_last_vertex(ctx->heat_path)->weight = T->value;
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-/*******************************************************************************
* Local functions
******************************************************************************/
res_T
@@ -135,7 +83,6 @@ XD(boundary_path)
}
if(res != RES_OK) goto error;
-#if 1
if(T->done) goto exit;
/* Handling limit boundary condition, i.e. the trajectory originates from a
@@ -155,11 +102,10 @@ XD(boundary_path)
* computationally, and in fact it's also a handy way of testing
* border cases */
if(T->func == XD(convective_path) || T->func == XD(conductive_path)) {
- res = XD(handle_known_medium_temperature)(scn, ctx, rwalk, T);
+ res = XD(query_medium_temperature_from_boundary)(scn, ctx, rwalk, T);
if(res != RES_OK) goto error;
if(T->done) goto exit; /* That's all folks */
}
-#endif
exit:
return res;
diff --git a/src/sdis_heat_path_boundary_Xd_c.h b/src/sdis_heat_path_boundary_Xd_c.h
@@ -1138,4 +1138,53 @@ XD(check_Tref)
return RES_OK;
}
+/* This function checks whether the random walk is on a boundary and, if so,
+ * verifies that the temperature of the medium attached to the interface is
+ * known. This medium can be different from the medium of the enclosure. Indeed,
+ * the enclosure can contain several media used to set the temperatures of
+ * several boundary conditions. Hence this function, which queries the medium on
+ * the path coming from a boundary */
+res_T
+XD(query_medium_temperature_from_boundary)
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk* rwalk,
+ struct temperature* T)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* mdm = NULL;
+ double temperature = SDIS_TEMPERATURE_NONE;
+ res_T res = RES_OK;
+ ASSERT(scn && ctx && rwalk && T);
+
+ /* Not at an interface */
+ if(SXD_HIT_NONE(&rwalk->XD(hit))) return RES_OK; /* Nothing to do */
+
+ interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id);
+ mdm = rwalk->hit_side==SDIS_FRONT ? interf->medium_front: interf->medium_back;
+
+ temperature = medium_get_temperature(mdm, &rwalk->vtx);
+
+ /* Check if the temperature is known */
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(temperature)) goto exit;
+
+ T->value += temperature;
+ T->done = 1;
+
+ if(ctx->green_path) {
+ res = green_path_set_limit_vertex
+ (ctx->green_path, mdm, &rwalk->vtx, rwalk->elapsed_time);
+ if(res != RES_OK) goto error;
+ }
+
+ if(ctx->heat_path) {
+ heat_path_get_last_vertex(ctx->heat_path)->weight = T->value;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
#include "sdis_Xd_end.h"
diff --git a/src/sdis_heat_path_boundary_Xd_handle_external_net_flux.h b/src/sdis_heat_path_boundary_Xd_handle_external_net_flux.h
@@ -503,6 +503,13 @@ XD(handle_external_net_flux)
handle_flux = handle_flux && (scn->source != NULL);
if(!handle_flux) goto exit;
+ /* Emissivity is null <=> external flux is null. Nothing to do */
+ src_id = sdis_source_get_id(scn->source);
+ emissivity = interface_side_get_emissivity(args->interf, src_id, &frag);
+ res = interface_side_check_emissivity(scn->dev, emissivity, frag.P, frag.time);
+ if(res != RES_OK) goto error;
+ if(emissivity == 0) goto exit;
+
/* Sample the external source */
res = source_sample
(scn->source, rng, frag.P, frag.time, &src_sample);
@@ -536,10 +543,6 @@ XD(handle_external_net_flux)
incident_flux_direct + incident_flux_diffuse.reflected;
/* Calculate the net flux [W/m^2] */
- src_id = sdis_source_get_id(scn->source);
- emissivity = interface_side_get_emissivity(args->interf, src_id, &frag);
- res = interface_side_check_emissivity(scn->dev, emissivity, frag.P, frag.time);
- if(res != RES_OK) goto error;
net_flux = incident_flux * emissivity; /* [W/m^2] */
/* Calculate the net flux from the radiance source scattered at least once by
diff --git a/src/sdis_heat_path_boundary_c.h b/src/sdis_heat_path_boundary_c.h
@@ -204,6 +204,24 @@ check_Tref_3d
const double Tref,
const char* call_func_name);
+/* Query medium temperature from the boundary. This medium can be different
+ * from the medium of the enclosure. Hence this function, which queries the
+ * medium on the path coming from a boundary. If the temperature is known, T is
+ * set to done. */
+extern LOCAL_SYM res_T
+query_medium_temperature_from_boundary_2d
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk* rwalk,
+ struct temperature* T);
+
+extern LOCAL_SYM res_T
+query_medium_temperature_from_boundary_3d
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk* rwalk,
+ struct temperature* T);
+
/*******************************************************************************
* Boundary sub-paths
******************************************************************************/
diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h
@@ -28,6 +28,84 @@
#include "sdis_Xd_begin.h"
/*******************************************************************************
+ * Non generic local functions
+ ******************************************************************************/
+#ifndef SDIS_HEAT_PATH_RADIATIVE_XD_H
+#define SDIS_HEAT_PATH_RADIATIVE_XD_H
+
+static res_T
+set_limit_radiative_temperature
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk* rwalk,
+ /* Direction along which the random walk reached the radiative environment */
+ const float dir[3],
+ const int branch_id,
+ struct temperature* T)
+{
+ struct sdis_radiative_ray ray = SDIS_RADIATIVE_RAY_NULL;
+ double trad = 0; /* [K] */
+ res_T res = RES_OK;
+
+ /* Check pre-conditions */
+ ASSERT(scn && ctx && rwalk && dir && T);
+ ASSERT(SXD_HIT_NONE(&rwalk->XD(hit)));
+
+ rwalk->hit_side = SDIS_SIDE_NULL__;
+ d3_set_f3(rwalk->dir, dir);
+ d3_normalize(rwalk->dir, rwalk->dir);
+ d3_set(ray.dir, rwalk->dir);
+
+ trad = radiative_env_get_temperature(scn->radenv, &ray);
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(trad)) {
+ log_err(scn->dev,
+ "%s:%s: the random walk has reached an invalid radiative environment from "
+ "position (%g, %g, %g) along direction (%g, %g, %g): the temperature is "
+ "unknown. This may be due to numerical inaccuracies or inconsistencies "
+ "in the simulated system (e.g. non-closed geometry). For systems where "
+ "sampled paths can reach such a temperature, we need to define a valid "
+ "radiative temperature, i.e. one with a known temperature.\n",
+ __FILE__, FUNC_NAME, SPLIT3(rwalk->vtx.P), SPLIT3(rwalk->dir));
+ res = RES_BAD_OP;
+ goto error;
+ }
+
+ /* The limit condition is reached */
+ T->value += trad;
+ T->done = 1;
+
+ /* Update green path */
+ if(ctx->green_path) {
+ res = green_path_set_limit_radiative_ray
+ (ctx->green_path, &ray, rwalk->elapsed_time);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Record the limit vertex of the sampled path. Set it arbitrarily at a
+ * distance of 0.1 meters from the surface along the direction reaching the
+ * radiative environment */
+ if(ctx->heat_path) {
+ const float empirical_dst = 0.1f * (float)scn->fp_to_meter;
+ struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
+
+ vtx = rwalk->vtx;
+ vtx.P[0] += dir[0] * empirical_dst;
+ vtx.P[1] += dir[1] * empirical_dst;
+ vtx.P[2] += dir[2] * empirical_dst;
+ res = register_heat_vertex(ctx->heat_path, &vtx, T->value,
+ SDIS_HEAT_VERTEX_RADIATIVE, branch_id);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#endif /* SDIS_HEAT_PATH_RADIATIVE_XD_H */
+
+/*******************************************************************************
* Local functions
******************************************************************************/
res_T
@@ -72,6 +150,8 @@ XD(trace_radiative_path)
/* Trace the radiative ray */
filter_data.XD(hit) = rwalk->XD(hit);
filter_data.epsilon = 1.e-6;
+ filter_data.scn = scn; /* Enable the filtering wrt the enclosure id */
+ filter_data.enc_id = rwalk->enc_id;
#if (SDIS_XD_DIMENSION == 2)
SXD(scene_view_trace_ray_3d
(scn->sXd(view), pos, dir, range, &filter_data, &rwalk->XD(hit)));
@@ -79,53 +159,14 @@ XD(trace_radiative_path)
SXD(scene_view_trace_ray
(scn->sXd(view), pos, dir, range, &filter_data, &rwalk->XD(hit)));
#endif
- if(SXD_HIT_NONE(&rwalk->XD(hit))) { /* Fetch the ambient radiative temperature */
- struct sdis_radiative_ray ray = SDIS_RADIATIVE_RAY_NULL;
- double trad = 0; /* [K] */
-
- rwalk->hit_side = SDIS_SIDE_NULL__;
- d3_set_f3(rwalk->dir, dir);
- d3_normalize(rwalk->dir, rwalk->dir);
- d3_set(ray.dir, rwalk->dir);
-
- trad = radiative_env_get_temperature(scn->radenv, &ray);
- if(SDIS_TEMPERATURE_IS_UNKNOWN(trad)) {
- log_err(scn->dev,
- "%s: the random walk has reached an invalid radiative environment from "
- "position `%g %g %g' along direction `%g %g %g': the temperature is "
- "unknown. This may be due to numerical inaccuracies or inconsistencies "
- "in the simulated system (e.g. non-closed geometry). For systems where "
- "random walks can reach such a temperature, we need to define a valid "
- "radiative temperature, i.e. one with a known temperature.\n",
- FUNC_NAME, SPLIT3(rwalk->vtx.P), SPLIT3(rwalk->dir));
- res = RES_BAD_OP;
- goto error;
- }
-
- T->value += trad;
- T->done = 1;
-
- if(ctx->green_path) {
- res = green_path_set_limit_radiative_ray
- (ctx->green_path, &ray, rwalk->elapsed_time);
- if(res != RES_OK) goto error;
- }
-
- if(ctx->heat_path) {
- const float empirical_dst = 0.1f * (float)scn->fp_to_meter;
- struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
-
- vtx = rwalk->vtx;
- vtx.P[0] += dir[0] * empirical_dst;
- vtx.P[1] += dir[1] * empirical_dst;
- vtx.P[2] += dir[2] * empirical_dst;
- res = register_heat_vertex(ctx->heat_path, &vtx, T->value,
- SDIS_HEAT_VERTEX_RADIATIVE, branch_id);
- if(res != RES_OK) goto error;
- }
-
- /* Stop the radiative path */
- break;
+
+ /* The path reaches the radiative environment */
+ if(SXD_HIT_NONE(&rwalk->XD(hit))) {
+ res = set_limit_radiative_temperature(scn, ctx, rwalk, dir, branch_id, T);
+ if(res != RES_OK) goto error;
+
+ ASSERT(T->done);
+ break; /* Stop the radiative path */
}
/* Define the hit side */
@@ -167,7 +208,12 @@ XD(trace_radiative_path)
fX(normalize)(N, rwalk->XD(hit).normal);
if(rwalk->hit_side == SDIS_BACK) fX(minus)(N, N);
- /* Check that the radiative path still lies in the same enclosure */
+ /* Check that the radiative path is still within the same enclosure. Note
+ * that this may not be the case, even if the filtering of intersections
+ * relative to the current enclosure is enabled. This filtering is only
+ * performed for intersections on a boundary between primitives. As a
+ * consequence, a threshold effect on how "intersections on a boundary" are
+ * detected could lead to this situation */
scene_get_enclosure_ids(scn, rwalk->XD(hit).prim.prim_id, enc_ids);
chk_enc_id = rwalk->hit_side == SDIS_FRONT ? enc_ids[0] : enc_ids[1];
if(chk_enc_id != rwalk->enc_id) {
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -15,6 +15,7 @@
#include "sdis_device_c.h"
#include "sdis_heat_path.h"
+#include "sdis_heat_path_boundary_c.h"
#include "sdis_interface_c.h"
#include "sdis_log.h"
#include "sdis_medium_c.h"
@@ -449,10 +450,28 @@ XD(boundary_flux_realisation)
/* Compute fluid temperature */
if(compute_convective) {
RESET_WALK(fluid_side, enc_ids[fluid_side]);
- T.func = XD(convective_path);
- res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T);
+
+ /* Check whether the temperature of the fluid is known by querying it from
+ * its boundary. This makes it possible to handle situations where fluids
+ * are used as Robin's boundary condition. In this case, a geometric
+ * enclosure may have several fluids, each defining the temperature of a
+ * boundary condition. Sampling of convective paths in such an enclosure is
+ * forbidden since this enclosure does not exist for this path space: it is
+ * beyond its boundary */
+ res = XD(query_medium_temperature_from_boundary)(scn, &ctx, &rwalk, &T);
if(res != RES_OK) return res;
- result->Tfluid = T.value;
+
+ /* Robin's boundary condition */
+ if(T.done) {
+ result->Tfluid = T.value;
+
+ /* Sample a convective path */
+ } else {
+ T.func = XD(convective_path);
+ res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T);
+ if(res != RES_OK) return res;
+ result->Tfluid = T.value;
+ }
}
#undef SET_PARAM
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -516,7 +516,24 @@ XD(hit_filter_function)
reject_hit = hit_shared_edge
(&hit_from->prim, &hit->prim, hit_from->uv, hit->uv, org, pos);
#endif
- return reject_hit;
+ if(reject_hit) return 1;
+ }
+
+ /* If the hit to be considered is (approximately) on a boundary between 2
+ * primitives, it may belong to the wrong primitive, i.e. the one that doesn't
+ * "face" the direction of the ray. We therefore check the enclosure towards
+ * which it is directed, and reject it if it is not the same as the one from
+ * which the ray originates */
+ if(filter_data->scn && HIT_ON_BOUNDARY(hit, org, dir)) {
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
+ unsigned chk_enc_id = ENCLOSURE_ID_NULL;
+
+ scene_get_enclosure_ids(filter_data->scn, hit->prim.prim_id, enc_ids);
+ chk_enc_id = fX(dot)(dir, hit->normal) < 0
+ ? enc_ids[0] /* Front */
+ : enc_ids[1]; /* Back */
+
+ if(chk_enc_id != filter_data->enc_id) return 1;
}
return 0;
}
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -40,6 +40,11 @@ struct hit_filter_data {
struct s3d_hit hit_3d;
double epsilon; /* Threshold defining roughly equal intersections */
+ /* When a scene is defined, primitives that do not point to the defined
+ * enclosure are filtered out */
+ struct sdis_scene* scn; /* NULL <=> do not filter wrt enc_id */
+ unsigned enc_id;
+
/* Bypass the regular filter function */
s2d_hit_filter_function_T custom_filter_2d;
s3d_hit_filter_function_T custom_filter_3d;
@@ -47,7 +52,8 @@ struct hit_filter_data {
/* Custom filter query data. It is ignored if custom_filter is NULL */
void* custom_filter_data;
};
-#define HIT_FILTER_DATA_NULL__ {S2D_HIT_NULL__,S3D_HIT_NULL__,0,NULL,NULL,NULL}
+#define HIT_FILTER_DATA_NULL__ \
+ {S2D_HIT_NULL__,S3D_HIT_NULL__,0,NULL,ENCLOSURE_ID_NULL,NULL,NULL,NULL}
static const struct hit_filter_data HIT_FILTER_DATA_NULL =
HIT_FILTER_DATA_NULL__;
diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c
@@ -266,7 +266,8 @@ main(int argc, char** argv)
{
struct sdis_data* data = NULL;
struct sdis_device* dev = NULL;
- struct sdis_medium* fluid = NULL;
+ struct sdis_medium* fluid1 = NULL;
+ struct sdis_medium* fluid2 = NULL;
struct sdis_medium* solid = NULL;
struct sdis_interface* interf_adiabatic = NULL;
struct sdis_interface* interf_Tb = NULL;
@@ -299,13 +300,19 @@ main(int argc, char** argv)
create_default_device(&argc, &argv, &is_master_process, &dev);
radenv = create_radenv(dev);
- /* Create the fluid medium */
+ /* Create the fluid medium. In fact, create two fluid media, even if they are
+ * identical, to check that Robin's boundary conditions can be defined with
+ * several media, without compromising path sampling. Convective paths cannot
+ * be sampled in enclosures with several media, but radiative paths can be,
+ * and it should be possible to calculate the boundary flux without any
+ * problem */
OK(sdis_data_create
(dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
fluid_args = sdis_data_get(data);
fluid_args->temperature = Tf;
fluid_shader.temperature = fluid_get_temperature;
- OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid));
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1));
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2));
OK(sdis_data_ref_put(data));
/* Create the solid_medium */
@@ -329,7 +336,7 @@ main(int argc, char** argv)
interf_props->temperature = SDIS_TEMPERATURE_NONE;
interf_props->emissivity = 0;
OK(sdis_interface_create
- (dev, solid, fluid, &interf_shader, data, &interf_adiabatic));
+ (dev, solid, fluid1, &interf_shader, data, &interf_adiabatic));
OK(sdis_data_ref_put(data));
/* Create the Tb interface */
@@ -342,7 +349,7 @@ main(int argc, char** argv)
interf_shader.back.emissivity = interface_get_emissivity;
interf_shader.back.reference_temperature = interface_get_reference_temperature;
OK(sdis_interface_create
- (dev, solid, fluid, &interf_shader, data, &interf_Tb));
+ (dev, solid, fluid1, &interf_shader, data, &interf_Tb));
interf_shader.back.emissivity = NULL;
OK(sdis_data_ref_put(data));
@@ -356,13 +363,14 @@ main(int argc, char** argv)
interf_shader.back.emissivity = interface_get_emissivity;
interf_shader.back.reference_temperature = interface_get_reference_temperature;
OK(sdis_interface_create
- (dev, solid, fluid, &interf_shader, data, &interf_H));
+ (dev, solid, fluid2, &interf_shader, data, &interf_H));
interf_shader.back.emissivity = NULL;
OK(sdis_data_ref_put(data));
/* Release the media */
OK(sdis_medium_ref_put(solid));
- OK(sdis_medium_ref_put(fluid));
+ OK(sdis_medium_ref_put(fluid1));
+ OK(sdis_medium_ref_put(fluid2));
/* Map the interfaces to their box triangles */
box_interfaces[0] = box_interfaces[1] = interf_adiabatic; /* Front */
diff --git a/src/test_sdis_unsteady_analytic_profile.c b/src/test_sdis_unsteady_analytic_profile.c
@@ -22,7 +22,7 @@
/*
* The system is an unsteady-state temperature profile, meaning that at any
- * point, at any: time, we can analytically calculate the temperature. We
+ * point, at any time, we can analytically calculate the temperature. We
* immerse in this temperature field a supershape representing a solid in which
* we want to evaluate the temperature by Monte Carlo at a given position and
* observation time. On the Monte Carlo side, the temperature of the supershape
