commit dd0b8a5f28901a711bc02284d5afc439f91e68f2
parent 37585c58e000230eef6c2462884240e100332286
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 16 Jan 2024 15:31:50 +0100
Add and test external flux support for green functions
In fact, green functions could be evaluated on scenes with an external
source, but this was simply ignored. This commit updates the external
flux calculation to allow it to be registered in the green function,
i.e. the external flux is calculated in relation to the source power.
This external flux term is registered in the green function, then
multiplied by the source power to obtain its contribution. The external
flux term is stored per path of the green function to facilitate
calculation of the standard error. Note that if the source's power is an
input parameter to the green function, the source's positions must
remain unchanged between its estimation and its evaluations. But the
solver can't verify this constraint, since the various source positions
are returned by a functor; to ensure that they remain unchanged, the
caller must add a signature to the green function based on the source
positions used to estimate the green function.
The API function sdis_green_path_get_external_flux_term has been added
to allow the caller to retrieve the external flux contribution by path.
In particular, it is used in tests to manually solve the green function
and compare the result obtained with the estimate returned by the
solver.
The external flux test has finally been updated to validate external
flux calculations using the green function: normal resolution is
performed for 2D scenes, while 3D scenes use green estimation. This
validates both 2D and 3D scenes, and the external flux calculation using
both the normal method and the green function, without prohibitively
increasing test time.
Diffstat:
10 files changed, 202 insertions(+), 38 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -1374,6 +1374,14 @@ sdis_green_path_for_each_flux_term
sdis_process_interface_flux_term_T func,
void* context);
+/* Return the external flux term, i.e. the relative net flux along the path from
+ * the external source. Multiply it by the power of the source to obtain its
+ * contribution to the path. */
+SDIS_API res_T
+sdis_green_path_get_external_flux_term
+ (struct sdis_green_path* path,
+ double* external_flux_term); /* [W/m^2] */
+
/*******************************************************************************
* Heat path API
******************************************************************************/
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -21,6 +21,7 @@
#include "sdis_medium_c.h"
#include "sdis_misc.h"
#include "sdis_scene_c.h"
+#include "sdis_source_c.h"
#include <star/ssp.h>
@@ -83,6 +84,7 @@ flux_term_init(struct mem_allocator* allocator, struct flux_term* term)
struct green_path {
double elapsed_time;
+ double external_flux_term; /* [W/m^2] */
struct darray_flux_term flux_terms; /* List of flux terms */
struct darray_power_term power_terms; /* List of volumic power terms */
union {
@@ -105,6 +107,7 @@ green_path_init(struct mem_allocator* allocator, struct green_path* path)
darray_flux_term_init(allocator, &path->flux_terms);
darray_power_term_init(allocator, &path->power_terms);
path->elapsed_time = -INF;
+ path->external_flux_term = 0;
path->limit.vertex = SDIS_RWALK_VERTEX_NULL;
path->limit.fragment = SDIS_INTERFACE_FRAGMENT_NULL;
path->limit_id = UINT_MAX;
@@ -127,6 +130,7 @@ green_path_copy(struct green_path* dst, const struct green_path* src)
res_T res = RES_OK;
ASSERT(dst && src);
dst->elapsed_time = src->elapsed_time;
+ dst->external_flux_term = src->external_flux_term;
dst->limit = src->limit;
dst->limit_id = src->limit_id;
dst->end_type = src->end_type;
@@ -145,6 +149,7 @@ green_path_copy_and_clear(struct green_path* dst, struct green_path* src)
res_T res = RES_OK;
ASSERT(dst && src);
dst->elapsed_time = src->elapsed_time;
+ dst->external_flux_term = src->external_flux_term;
dst->limit = src->limit;
dst->limit_id = src->limit_id;
dst->end_type = src->end_type;
@@ -164,6 +169,7 @@ green_path_copy_and_release(struct green_path* dst, struct green_path* src)
res_T res = RES_OK;
ASSERT(dst && src);
dst->elapsed_time = src->elapsed_time;
+ dst->external_flux_term = src->external_flux_term;
dst->limit = src->limit;
dst->limit_id = src->limit_id;
dst->end_type = src->end_type;
@@ -192,6 +198,7 @@ green_path_write(const struct green_path* path, FILE* stream)
/* Write elapsed time */
WRITE(&path->elapsed_time, 1);
+ WRITE(&path->external_flux_term, 1);
/* Write the list of flux terms */
sz = darray_flux_term_size_get(&path->flux_terms);
@@ -242,6 +249,7 @@ green_path_read(struct green_path* path, FILE* stream)
/* Read elapsed time */
READ(&path->elapsed_time, 1);
+ READ(&path->external_flux_term, 1);
/* Read the list of flux terms */
READ(&sz, 1);
@@ -397,7 +405,7 @@ green_function_solve_path
const size_t ipath,
double* weight)
{
- struct sdis_ambient_radiative_temperature trad =
+ struct sdis_ambient_radiative_temperature trad =
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
const struct power_term* power_terms = NULL;
const struct flux_term* flux_terms = NULL;
@@ -409,6 +417,7 @@ green_function_solve_path
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
double power;
double flux;
+ double external_flux;
double end_temperature;
size_t i, n;
res_T res = RES_OK;
@@ -420,7 +429,7 @@ green_function_solve_path
goto error;
}
- /* Compute medium power terms */
+ /* Compute medium powers */
power = 0;
n = darray_power_term_size_get(&path->power_terms);
power_terms = darray_power_term_cdata_get(&path->power_terms);
@@ -441,6 +450,13 @@ green_function_solve_path
flux += flux_terms[i].term * interface_side_get_flux(interf, &frag);
}
+ /* Compute external flux */
+ external_flux = 0;
+ if(green->scn->source) {
+ external_flux =
+ path->external_flux_term * source_get_power(green->scn->source);
+ }
+
/* Compute path's end temperature */
switch(path->end_type) {
case SDIS_GREEN_PATH_END_AT_INTERFACE:
@@ -466,7 +482,7 @@ green_function_solve_path
}
/* Compute the path weight */
- *weight = power + flux + end_temperature;
+ *weight = power + flux + external_flux + end_temperature;
exit:
return res;
@@ -1032,7 +1048,7 @@ error:
res_T
sdis_green_function_get_scene
- (const struct sdis_green_function* green,
+ (const struct sdis_green_function* green,
struct sdis_scene** scn)
{
if(!green || !scn) return RES_BAD_ARG;
@@ -1352,6 +1368,34 @@ error:
goto exit;
}
+res_T
+sdis_green_path_get_external_flux_term
+ (struct sdis_green_path* path_handle,
+ double* external_flux_term)
+{
+ const struct green_path* path = NULL;
+ struct sdis_green_function* green = NULL;
+ res_T res = RES_OK;
+
+ if(!path_handle || !external_flux_term) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ green = path_handle->green__; (void)green;
+ ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths));
+
+ path = darray_green_path_cdata_get(&green->paths) + path_handle->id__;
+
+ *external_flux_term = path->external_flux_term;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -1710,3 +1754,13 @@ exit:
error:
goto exit;
}
+
+res_T
+green_path_add_external_flux_term
+ (struct green_path_handle* handle,
+ const double val) /* [W/m^2/sr] */
+{
+ ASSERT(handle);
+ handle->path->external_flux_term += val;
+ return RES_OK;
+}
diff --git a/src/sdis_green.h b/src/sdis_green.h
@@ -20,9 +20,9 @@
#include <rsys/rsys.h>
/* Current version the green function data structure. One should increment it
- * and perform a version management onto serialized data when the gren function
+ * and perform a version management onto serialized data when the green function
* data structure is updated. */
-static const int SDIS_GREEN_FUNCTION_VERSION = 2;
+static const int SDIS_GREEN_FUNCTION_VERSION = 3;
/* Forward declaration */
struct accum;
@@ -106,5 +106,10 @@ green_path_add_flux_term
const struct sdis_interface_fragment* fragment,
const double term);
+extern LOCAL_SYM res_T
+green_path_add_external_flux_term
+ (struct green_path_handle* handle,
+ const double term); /* [W/m^2/sr] */
+
#endif /* SDIS_GREEN_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
@@ -240,7 +240,12 @@ XD(direct_contribution)
XD(trace_ray)(scn, pos, sample->dir, sample->dst, hit_from, &hit);
if(!SXD_HIT_NONE(&hit)) return 0; /* [W/m^2/sr] */
- return sample->radiance; /* [W/m^2/sr] */
+ /* Note that the value returned is not the source's actual radiance, but the
+ * radiance relative to the source's power. Care must therefore be taken to
+ * multiply it by the power of the source to obtain its real contribution.
+ * This trick makes it possible to manage the external flux in the green
+ * function. */
+ return sample->radiance_term; /* [W/m^2/sr] */
}
static INLINE void
@@ -430,6 +435,7 @@ XD(handle_external_net_flux)
double incident_flux_diffuse = 0; /* [W/m^2] */
double incident_flux_direct = 0; /* [W/m^2] */
double net_flux = 0; /* [W/m^2] */
+ double external_flux_term = 0; /* [W/m^2] */
/* Sampled path */
double N[3] = {0}; /* Normal. Always in 3D */
@@ -493,8 +499,19 @@ XD(handle_external_net_flux)
if(res != RES_OK) goto error;
net_flux = incident_flux * emissivity; /* [W/m^2] */
+ /* Until now, the net flux was calculated in relation to the source power.
+ * What is calculated is the external flux term of the green function. This
+ * must be multiplied by the source power to obtain the actual external flux*/
+ external_flux_term = net_flux / (args->h_radi + args->h_conv + args->h_cond);
+
/* Update the Monte Carlo weight */
- T->value += net_flux / (args->h_radi + args->h_conv + args->h_cond);
+ T->value += external_flux_term * source_get_power(scn->source);
+
+ /* Register the external net flux term */
+ if(args->green_path) {
+ res = green_path_add_external_flux_term(args->green_path, external_flux_term);
+ if(res != RES_OK) goto error;
+ }
exit:
return res;
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h
@@ -246,6 +246,7 @@ XD(solid_fluid_boundary_picard1_path)
handle_external_net_flux_args.interf = interf;
handle_external_net_flux_args.frag = frag;
handle_external_net_flux_args.hit = &rwalk->hit;
+ handle_external_net_flux_args.green_path = ctx->green_path;
handle_external_net_flux_args.picard_order = get_picard_order(ctx);
handle_external_net_flux_args.h_cond = h_cond;
handle_external_net_flux_args.h_conv = h_conv;
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -482,6 +482,13 @@ scene_compute_hash(const struct sdis_scene* scn, hash256_T hash)
SHA256_UPD(&scn->trad.reference, 1);
SHA256_UPD(&scn->tmax, 1);
SHA256_UPD(&scn->fp_to_meter, 1);
+
+ if(scn->source) {
+ hash256_T src_hash;
+ source_compute_signature(scn->source, src_hash);
+ sha256_ctx_update(&sha256_ctx, src_hash, sizeof(hash256_T));
+ }
+
FOR_EACH(iprim, 0, nprims) {
struct sdis_interface* interf = NULL;
size_t ivert;
diff --git a/src/sdis_source.c b/src/sdis_source.c
@@ -186,7 +186,9 @@ source_sample
d3_set(sample->dir, main_dir);
sample->pdf = 1;
sample->dst = dst;
- sample->radiance = src->spherical.power / (4*PI*dst*dst); /* [W/m^2/sr] */
+ sample->power = src->spherical.power; /* [W] */
+ sample->radiance_term = 1.0 / (4*PI*dst*dst); /* [W/m^2/sr] */
+ sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
/* Spherical source */
} else {
@@ -197,10 +199,11 @@ source_sample
ssp_ran_sphere_cap_uniform /* pdf = 1/(2*PI*(1-cos(half_angle))) */
(rng, main_dir, cos_half_angle, sample->dir, &sample->pdf);
-
/* Set other sample variables */
sample->dst = dst - radius; /* From pos to source boundaries [m] */
- sample->radiance = src->spherical.power / (PI*area); /* [W/m^2/sr] */
+ sample->power = src->spherical.power; /* [W] */
+ sample->radiance_term = 1.0 / (PI*area); /* [W/m^2/sr] */
+ sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
}
exit:
@@ -257,7 +260,7 @@ source_trace_to
goto error;
}
- /* Compute the the half angle of the source as seen from pos */
+ /* Compute the half angle of the source as seen from pos */
half_angle = asin(radius/dst);
/* The source is missed */
@@ -271,7 +274,9 @@ source_trace_to
d3_set(sample->dir, dir);
sample->pdf = 1;
sample->dst = dst - radius; /* From pos to source boundaries [m] */
- sample->radiance = src->spherical.power / (PI*area); /* [W/m^2/sr] */
+ sample->power = src->spherical.power; /* [W] */
+ sample->radiance_term = 1.0 / (PI*area); /* [W/m^2/sr] */
+ sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
}
exit:
@@ -287,3 +292,20 @@ source_get_power(const struct sdis_source* src)
ASSERT(src);
return src->spherical.power;
}
+
+void
+source_compute_signature(const struct sdis_source* src, hash256_T hash)
+{
+ struct sha256_ctx ctx;
+ ASSERT(src && hash);
+
+ /* Calculate the source signature. Currently, it is only the source radius.
+ * But the Source API is designed to be independent of source type. In the
+ * future, the source will not necessarily be spherical, so the data to be
+ * hashed will depend on the type of source. This function anticipate this by
+ * calculating a hash even if it is currently dispensable. */
+ sha256_ctx_init(&ctx);
+ sha256_ctx_update
+ (&ctx, (const char*)&src->spherical.radius, sizeof(src->spherical.radius));
+ sha256_ctx_finalize(&ctx, hash);
+}
diff --git a/src/sdis_source_c.h b/src/sdis_source_c.h
@@ -16,6 +16,7 @@
#ifndef SDIS_SOURCE_C_H
#define SDIS_SOURCE_C_H
+#include <rsys/hash.h>
#include <rsys/rsys.h>
struct sdis_source;
@@ -25,9 +26,17 @@ struct source_sample {
double dir[3]; /* Direction _to_ the source */
double pdf; /* pdf of sampled direction */
double dst; /* Distance to the source [m] */
+ double power; /* [W] */
double radiance; /* [W/m^2/sr] */
+
+ /* Radiance relative to power, i.e. the source power is assumed to be equal to
+ * 1. It must be multiplied by the source power to obtain the actual radiance
+ * of the source. In other words, this variable defines the contribution of
+ * the source independently of its power, and can therefore be recorded in the
+ * green function */
+ double radiance_term; /* [W/m^2/sr] */
};
-#define SOURCE_SAMPLE_NULL__ {{0,0,0}, 0, 0, 0}
+#define SOURCE_SAMPLE_NULL__ {{0,0,0}, 0, 0, 0, 0, 0}
static const struct source_sample SOURCE_SAMPLE_NULL = SOURCE_SAMPLE_NULL__;
/* Helper macro used to define whether a sample is valid or not */
@@ -56,4 +65,9 @@ extern LOCAL_SYM double /* [W] */
source_get_power
(const struct sdis_source* source);
+extern LOCAL_SYM void
+source_compute_signature
+ (const struct sdis_source* source,
+ hash256_T hash);
+
#endif /* SDIS_SOURCE_C_H */
diff --git a/src/test_sdis_external_flux.c b/src/test_sdis_external_flux.c
@@ -422,8 +422,6 @@ check
struct sdis_solve_probe_args probe_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
struct sdis_mc T = SDIS_MC_NULL;
struct sdis_estimator* estimator = NULL;
- enum sdis_scene_dimension dim;
- const char* dim_str = NULL;
probe_args.position[0] = -0.05;
probe_args.position[1] = 1000;
@@ -432,14 +430,37 @@ check
OK(sdis_solve_probe(scn, &probe_args, &estimator));
OK(sdis_estimator_get_temperature(estimator, &T));
- OK(sdis_scene_get_dimension(scn, &dim));
- switch(dim) {
- case SDIS_SCENE_2D: dim_str = "2D"; break;
- case SDIS_SCENE_3D: dim_str = "3D"; break;
- default: FATAL("Unreachable code.\n"); break;
- }
- printf("%s: T(%g, %g, %g) = %g ~ %g +/- %g\n",
- dim_str, SPLIT3(probe_args.position), analytical_ref, T.E, T.SE);
+ printf("T(%g, %g, %g) = %g ~ %g +/- %g\n",
+ SPLIT3(probe_args.position), analytical_ref, T.E, T.SE);
+ OK(sdis_estimator_ref_put(estimator));
+
+ CHK(eq_eps(analytical_ref, T.E, 3*T.SE));
+}
+
+static void
+check_green
+ (struct sdis_scene* scn,
+ const size_t nrealisations,
+ const double analytical_ref)
+{
+ struct sdis_solve_probe_args probe_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_green_function* green = NULL;
+ struct sdis_estimator* estimator = NULL;
+
+ probe_args.position[0] = -0.05;
+ probe_args.position[1] = 1000;
+ probe_args.position[2] = 0;
+ probe_args.nrealisations = nrealisations;
+ OK(sdis_solve_probe_green_function(scn, &probe_args, &green));
+ OK(sdis_green_function_solve(green, &estimator));
+ check_green_function(green);
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+
+ printf("T(%g, %g, %g) = %g ~ %g +/- %g\n",
+ SPLIT3(probe_args.position), analytical_ref, T.E, T.SE);
+ OK(sdis_green_function_ref_put(green));
OK(sdis_estimator_ref_put(estimator));
CHK(eq_eps(analytical_ref, T.E, 3*T.SE));
@@ -477,11 +498,11 @@ main(int argc, char** argv)
ground_interf_data->specular_fraction = 0; /* Lambertian */
check(scn_2d, 10000/* #réalisations */, 375.88/* Reference [K] */);
- check(scn_3d, 10000/* #réalisations */, 375.88/* Reference [K] */);
+ check_green(scn_3d, 10000/* #réalisations */, 375.88/* Reference [K] */);
ground_interf_data->specular_fraction = 1; /* Specular */
check(scn_2d, 100000/* #réalisations */, 417.77/* Reference [K] */);
- check(scn_3d, 100000/* #réalisations */, 417.77/* Reference [K] */);
+ check_green(scn_3d, 100000/* #réalisations */, 417.77/* Reference [K] */);
OK(sdis_device_ref_put(dev));
OK(sdis_medium_ref_put(fluid));
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -92,18 +92,33 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
struct sdis_solid_shader solid = SDIS_SOLID_SHADER_NULL;
struct sdis_fluid_shader fluid = SDIS_FLUID_SHADER_NULL;
struct sdis_interface_shader interf = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_green_function* green = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_source* source = NULL;
struct green_accum* acc = NULL;
struct sdis_data* data = NULL;
enum sdis_medium_type type;
enum sdis_green_path_end_type end_type;
+ double source_power = 0; /* [W] */
double power = 0;
double flux = 0;
+ double external_flux = 0; /* [W/m^2] */
double time, temp = 0;
double weight = 0;
CHK(path && ctx);
acc = ctx;
+ BA(sdis_green_path_get_green_function(NULL, NULL));
+ BA(sdis_green_path_get_green_function(path, NULL));
+ BA(sdis_green_path_get_green_function(NULL, &green));
+ OK(sdis_green_path_get_green_function(path, &green));
+
+ BA(sdis_green_function_get_scene(NULL, NULL));
+ BA(sdis_green_function_get_scene(NULL, &scn));
+ BA(sdis_green_function_get_scene(green, NULL));
+ OK(sdis_green_function_get_scene(green, &scn));
+
BA(sdis_green_path_for_each_power_term(NULL, accum_power_terms, &power));
BA(sdis_green_path_for_each_power_term(path, NULL, &acc));
OK(sdis_green_path_for_each_power_term(path, accum_power_terms, &power));
@@ -117,6 +132,17 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
BA(sdis_green_path_get_elapsed_time(NULL, &time));
OK(sdis_green_path_get_elapsed_time(path, &time));
+ BA(sdis_green_path_get_external_flux_term(NULL, &external_flux));
+ BA(sdis_green_path_get_external_flux_term(path, NULL));
+ OK(sdis_green_path_get_external_flux_term(path, &external_flux));
+ OK(sdis_scene_get_source(scn, &source));
+ if(source == NULL) {
+ CHK(external_flux == 0);
+ } else {
+ OK(sdis_source_get_power(source, &source_power));
+ external_flux *= source_power;
+ }
+
BA(sdis_green_path_get_end_type(NULL, NULL));
BA(sdis_green_path_get_end_type(path, NULL));
BA(sdis_green_path_get_end_type(NULL, &end_type));
@@ -127,18 +153,7 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
BA(sdis_green_path_get_limit_point(path, NULL));
if(end_type == SDIS_GREEN_PATH_END_RADIATIVE) {
struct sdis_ambient_radiative_temperature trad;
- struct sdis_green_function* green;
- struct sdis_scene* scn;
BO(sdis_green_path_get_limit_point(path, &pt));
- BA(sdis_green_path_get_green_function(NULL, NULL));
- BA(sdis_green_path_get_green_function(path, NULL));
- BA(sdis_green_path_get_green_function(NULL, &green));
- OK(sdis_green_path_get_green_function(path, &green));
-
- BA(sdis_green_function_get_scene(NULL, NULL));
- BA(sdis_green_function_get_scene(NULL, &scn));
- BA(sdis_green_function_get_scene(green, NULL));
- OK(sdis_green_function_get_scene(green, &scn));
BA(sdis_scene_get_ambient_radiative_temperature(NULL, NULL));
BA(sdis_scene_get_ambient_radiative_temperature(scn, NULL));
@@ -172,7 +187,7 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
}
}
- weight = temp + power + flux;
+ weight = temp + power + external_flux + flux;
acc->sum += weight;
acc->sum2 += weight*weight;
