commit a59ae312ec435aac26f68d910a026f7e82ca9766
parent 35309912f04233443f9dbe897286cb22fc9c4df2
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 18 Apr 2024 12:05:20 +0200
Merge branch 'feature_extern_diffuse_radiance' into develop
Diffstat:
10 files changed, 871 insertions(+), 113 deletions(-)
diff --git a/Makefile b/Makefile
@@ -183,6 +183,7 @@ TEST_SRC =\
src/test_sdis_data.c\
src/test_sdis_draw_external_flux.c\
src/test_sdis_enclosure_limit_conditions.c\
+ src/test_sdis_external_flux_with_diffuse_radiance.c\
src/test_sdis_flux.c\
src/test_sdis_flux2.c\
src/test_sdis_flux_with_h.c\
@@ -392,6 +393,7 @@ test_sdis_volumic_power4 \
# Tests based on Star-3DUT
################################################################################
src/test_sdis_draw_external_flux.d \
+src/test_sdis_external_flux_with_diffuse_radiance.d \
src/test_sdis_solid_random_walk_robustness.d \
src/test_sdis_solve_probe3.d \
src/test_sdis_unsteady_analytic_profile.d \
@@ -399,6 +401,7 @@ src/test_sdis_unsteady_analytic_profile.d \
@$(CC) $(TEST_CFLAGS) $(S3DUT_CFLAGS) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
src/test_sdis_draw_external_flux.o \
+src/test_sdis_external_flux_with_diffuse_radiance.o \
src/test_sdis_solid_random_walk_robustness.o \
src/test_sdis_solve_probe3.o \
src/test_sdis_unsteady_analytic_profile.o \
@@ -406,6 +409,7 @@ src/test_sdis_unsteady_analytic_profile.o \
$(CC) $(TEST_CFLAGS) $(S3DUT_CFLAGS) -c $(@:.o=.c) -o $@
test_sdis_draw_external_flux \
+test_sdis_external_flux_with_diffuse_radiance \
test_sdis_solid_random_walk_robustness \
test_sdis_solve_probe3 \
test_sdis_unsteady_analytic_profile \
diff --git a/src/sdis.h b/src/sdis.h
@@ -171,24 +171,39 @@ static const struct sdis_info SDIS_INFO_NULL = SDIS_INFO_NULL__;
/* Type of functor used to retrieve the source's position relative to time */
typedef void
(*sdis_get_position_T)
- (const double time,
+ (const double time, /* [s] */
double pos[3],
struct sdis_data* data);
/* Type of functor used to retrieve the source's power relative to time */
typedef double
(*sdis_get_power_T)
- (const double time,
+ (const double time, /* [s] */
+ struct sdis_data* data);
+
+/* Type of functor used to retrieve the diffuse part of the external radiance */
+typedef double /* [W/perpendicular m^2/sr] */
+(*sdis_get_diffuse_radiance_T)
+ (const double time, /* [s] */
+ const double dir[3],
struct sdis_data* data);
/* Input arguments of the sdis_spherical_source_create function */
struct sdis_spherical_source_create_args {
- sdis_get_position_T position; /* [m] */
+ sdis_get_position_T position; /* [m/fp_to_meter] */
sdis_get_power_T power; /* Total power [W] */
+
+ /* Describes the diffuse part of the source's radiance, i.e. the radiance
+ * emitted by the source and scattered at least once in the environment. This
+ * parameter is actually used to approximate a semi-transparent medium. Its
+ * value can be NULL, meaning that the source has not been scattered by the
+ * environment, or, to put it another way, that the source is in a vacuum. */
+ sdis_get_diffuse_radiance_T diffuse_radiance; /* [W/m^2/sr] */
+
struct sdis_data* data; /* Data sent to the position functor */
double radius; /* [m] */
};
-#define SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL__ {NULL, NULL, 0, 0}
+#define SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL__ {NULL, NULL, NULL, 0, 0}
static const struct sdis_spherical_source_create_args
SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL =
SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL__;
@@ -441,27 +456,49 @@ struct sdis_green_path_end {
static const struct sdis_green_path_end SDIS_GREEN_PATH_END_NULL =
SDIS_GREEN_PATH_END_NULL__;
-/* Functor used to process the paths registered against the green function */
+struct sdis_green_external_flux_terms {
+ /* Term relative to source power [K/W] */
+ double term_wrt_power;
+
+ /* Term relative to diffuse source radiance [K/W/m^2/sr] */
+ double term_wrt_diffuse_radiance;
+
+ double time; /* [s] */
+ double dir[3]; /* Direction on which term_wrt_diffuse_radiance depends */
+};
+#define SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL__ {0,0,0,{0,0,0}}
+static const struct sdis_green_external_flux_terms
+SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL = SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL__;
+
+/* Function profile used to process the paths stored in the green function */
typedef res_T
(*sdis_process_green_path_T)
(struct sdis_green_path* path,
void* context);
-/* Functor used to process the power factor registered along a green path for a
- * given medium */
+/* Function profile used to process power factors registered along a green path
+ * for a given medium */
typedef res_T
(*sdis_process_medium_power_term_T)
(struct sdis_medium* medium,
- const double power_term,
+ const double power_term, /* [K/W] */
void* context);
-/* Functor used to process the flux factor registered along a green path for a
- * given interface side */
+/* Function profile used to process flux factors recorded along a green path for
+ * a given interface side */
typedef res_T
(*sdis_process_interface_flux_term_T)
(struct sdis_interface* interf,
const enum sdis_side side,
- const double flux_term,
+ const double flux_term, /* [K/W/m^2] */
+ void* context);
+
+/* Function profile used to process external flux factors recorded along a green
+ * path */
+typedef res_T
+(*sdis_process_external_flux_terms_T)
+ (struct sdis_source* source,
+ const struct sdis_green_external_flux_terms* terms,
void* context);
/*******************************************************************************
@@ -1118,11 +1155,18 @@ SDIS_API res_T
sdis_source_ref_put
(struct sdis_source* source);
-SDIS_API double
+SDIS_API double /* [W] */
sdis_source_get_power
(struct sdis_source* source,
const double time); /* [s] */
+/* Return the source radiance that is diffused in the environment */
+SDIS_API double /* [W/m^2/sr*] */
+sdis_source_get_diffuse_radiance
+ (struct sdis_source* source,
+ const double time, /* [s] */
+ const double dir[3]);
+
SDIS_API unsigned
sdis_source_get_id
(const struct sdis_source* source);
@@ -1486,6 +1530,11 @@ sdis_green_function_get_flux_terms_count
(const struct sdis_green_path* path,
size_t* nterms);
+SDIS_API res_T
+sdis_green_function_get_external_flux_terms_count
+ (const struct sdis_green_path* path,
+ size_t* nterms);
+
/* Iterate over all "power terms" associated to the path. Multiply each term
* by the power of their associated medium, that is assumed to be constant in
* time and space, gives the medium power registered along the path. */
@@ -1504,13 +1553,12 @@ 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. */
+/* Iterate over all external flux terms associated to the path */
SDIS_API res_T
-sdis_green_path_get_external_flux_term
+sdis_green_path_for_each_external_flux_terms
(struct sdis_green_path* path,
- double* external_flux_term); /* [W/m^2] */
+ sdis_process_external_flux_terms_T func,
+ void* context);
/*******************************************************************************
* Heat path API
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -26,6 +26,7 @@
#include <star/ssp.h>
+#include <rsys/cstr.h>
#include <rsys/dynamic_array.h>
#include <rsys/hash_table.h>
#include <rsys/mem_allocator.h>
@@ -83,9 +84,24 @@ flux_term_init(struct mem_allocator* allocator, struct flux_term* term)
#define DARRAY_FUNCTOR_INIT flux_term_init
#include <rsys/dynamic_array.h>
+static INLINE void
+extflux_terms_init
+ (struct mem_allocator* allocator,
+ struct sdis_green_external_flux_terms* terms)
+{
+ ASSERT(terms); (void)allocator;
+ *terms = SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL;
+}
+
+/* Generate the dynamic array of the external flux terms */
+#define DARRAY_NAME extflux_terms
+#define DARRAY_DATA struct sdis_green_external_flux_terms
+#define DARRAY_FUNCTOR_INIT extflux_terms_init
+#include <rsys/dynamic_array.h>
+
struct green_path {
double elapsed_time;
- double external_flux_term; /* [W/m^2] */
+ struct darray_extflux_terms extflux_terms; /* List of external flux terms */
struct darray_flux_term flux_terms; /* List of flux terms */
struct darray_power_term power_terms; /* List of volumic power terms */
union {
@@ -106,10 +122,10 @@ static INLINE void
green_path_init(struct mem_allocator* allocator, struct green_path* path)
{
ASSERT(path);
+ darray_extflux_terms_init(allocator, &path->extflux_terms);
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.ray = SDIS_RADIATIVE_RAY_NULL;
@@ -125,6 +141,7 @@ green_path_release(struct green_path* path)
ASSERT(path);
darray_flux_term_release(&path->flux_terms);
darray_power_term_release(&path->power_terms);
+ darray_extflux_terms_release(&path->extflux_terms);
}
static INLINE res_T
@@ -133,12 +150,13 @@ 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;
dst->ilast_medium = src->ilast_medium;
dst->ilast_interf = src->ilast_interf;
+ res = darray_extflux_terms_copy(&dst->extflux_terms, &src->extflux_terms);
+ if(res != RES_OK) return res;
res = darray_flux_term_copy(&dst->flux_terms, &src->flux_terms);
if(res != RES_OK) return res;
res = darray_power_term_copy(&dst->power_terms, &src->power_terms);
@@ -152,12 +170,14 @@ 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;
dst->ilast_medium = src->ilast_medium;
dst->ilast_interf = src->ilast_interf;
+ res = darray_extflux_terms_copy_and_clear
+ (&dst->extflux_terms, &src->extflux_terms);
+ if(res != RES_OK) return res;
res = darray_flux_term_copy_and_clear(&dst->flux_terms, &src->flux_terms);
if(res != RES_OK) return res;
res = darray_power_term_copy_and_clear(&dst->power_terms, &src->power_terms);
@@ -172,12 +192,14 @@ 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;
dst->ilast_medium = src->ilast_medium;
dst->ilast_interf = src->ilast_interf;
+ res = darray_extflux_terms_copy_and_release
+ (&dst->extflux_terms, &src->extflux_terms);
+ if(res != RES_OK) return res;
res = darray_flux_term_copy_and_release(&dst->flux_terms, &src->flux_terms);
if(res != RES_OK) return res;
res = darray_power_term_copy_and_release(&dst->power_terms, &src->power_terms);
@@ -201,7 +223,11 @@ 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 external flux terms */
+ sz = darray_extflux_terms_size_get(&path->extflux_terms);
+ WRITE(&sz, 1);
+ WRITE(darray_extflux_terms_cdata_get(&path->extflux_terms), sz);
/* Write the list of flux terms */
sz = darray_flux_term_size_get(&path->flux_terms);
@@ -252,7 +278,12 @@ 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 external flux terms */
+ READ(&sz, 1);
+ res = darray_extflux_terms_resize(&path->extflux_terms, sz);
+ if(res != RES_OK) goto error;
+ READ(darray_extflux_terms_data_get(&path->extflux_terms), sz);
/* Read the list of flux terms */
READ(&sz, 1);
@@ -402,14 +433,105 @@ green_function_fetch_interf
return *pinterf;
}
+static double /* [K] */
+green_path_power_contribution
+ (struct sdis_green_function* green,
+ const struct green_path* path)
+{
+ double temperature = 0; /* [K] */
+ size_t i=0, n=0;
+
+ ASSERT(green && path);
+
+ n = darray_power_term_size_get(&path->power_terms);
+ FOR_EACH(i, 0, n) {
+ struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
+ const struct power_term* power_term = NULL;
+ const struct sdis_medium* medium = NULL;
+ double power = 0; /* [W] */
+
+ power_term = darray_power_term_cdata_get(&path->power_terms) + i;
+ medium = green_function_fetch_medium(green, power_term->id);
+
+ /* Dummy argument used only to satisfy the function profile used to recover
+ * power. Its position is unused, since power is assumed to be constant in
+ * space, and its time is set to infinity, since the green function is
+ * assumed to be evaluated at steady state */
+ vtx.time = INF;
+ power = solid_get_volumic_power(medium, &vtx);
+
+ temperature += power_term->term * power; /* [K] */
+ }
+ return temperature; /* [K] */
+}
+
+static double /* [K] */
+green_path_flux_contribution
+ (struct sdis_green_function* green,
+ const struct green_path* path)
+{
+ double temperature = 0;
+ size_t i=0, n=0;
+
+ ASSERT(green && path);
+
+ n = darray_flux_term_size_get(&path->flux_terms);
+ FOR_EACH(i, 0, n) {
+ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+ const struct flux_term* flux_term = NULL;
+ const struct sdis_interface* interf = NULL;
+ double flux = 0; /* [W/m^2] */
+
+ flux_term = darray_flux_term_cdata_get(&path->flux_terms) + i;
+ interf = green_function_fetch_interf(green, flux_term->id);
+
+ /* Interface fragment. Its position is unused, since flux is assumed to be
+ * constant in space, and its time is set to infinity, since the green
+ * function is assumed to be evaluated at steady state */
+ frag.time = INF;
+ frag.side = flux_term->side;
+ flux = interface_side_get_flux(interf, &frag);
+
+ temperature += flux_term->term * flux; /* [K] */
+ }
+ return temperature; /* [K] */
+}
+
+static double /* [K] */
+green_path_external_flux_contribution
+ (struct sdis_green_function* green,
+ const struct green_path* path)
+{
+ const struct sdis_source* extsrc = NULL;
+ double value = 0;
+ size_t i=0, n=0;
+
+ ASSERT(green && path);
+
+ if((extsrc = green->scn->source) == NULL) return 0;
+
+ n = darray_extflux_terms_size_get(&path->extflux_terms);
+ FOR_EACH(i, 0, n) {
+ const struct sdis_green_external_flux_terms* extflux = NULL;
+ double power = 0; /* [W] */
+ double diffrad = 0; /* [W/m^2/sr] */
+
+ extflux = darray_extflux_terms_cdata_get(&path->extflux_terms)+i;
+ power = source_get_power(extsrc, extflux->time);
+ diffrad = source_get_diffuse_radiance(extsrc, extflux->time, extflux->dir);
+
+ value += extflux->term_wrt_power * power; /* [K] */
+ value += extflux->term_wrt_diffuse_radiance * diffrad; /* [K] */
+ }
+ return value; /* [K] */
+}
+
static res_T
green_function_solve_path
(struct sdis_green_function* green,
const size_t ipath,
double* weight)
{
- const struct power_term* power_terms = NULL;
- const struct flux_term* flux_terms = NULL;
const struct green_path* path = NULL;
const struct sdis_medium* medium = NULL;
const struct sdis_interface* interf = NULL;
@@ -421,7 +543,6 @@ green_function_solve_path
double flux;
double external_flux;
double end_temperature;
- size_t i, n;
res_T res = RES_OK;
ASSERT(green && ipath < darray_green_path_size_get(&green->paths) && weight);
@@ -431,35 +552,9 @@ green_function_solve_path
goto error;
}
- /* Compute medium powers */
- power = 0;
- n = darray_power_term_size_get(&path->power_terms);
- power_terms = darray_power_term_cdata_get(&path->power_terms);
- FOR_EACH(i, 0, n) {
- vtx.time = INF;
- medium = green_function_fetch_medium(green, power_terms[i].id);
- power += power_terms[i].term * solid_get_volumic_power(medium, &vtx);
- }
-
- /* Compute interface fluxes */
- flux = 0;
- n = darray_flux_term_size_get(&path->flux_terms);
- flux_terms = darray_flux_term_cdata_get(&path->flux_terms);
- FOR_EACH(i, 0, n) {
- frag.time = INF;
- frag.side = flux_terms[i].side;
- interf = green_function_fetch_interf(green, flux_terms[i].id);
- flux += flux_terms[i].term * interface_side_get_flux(interf, &frag);
- }
-
- /* Compute external flux */
- external_flux = 0;
- if(green->scn->source) {
- /* NOTE: The power of the source is assumed to be constant over time and is
- * therefore recovered at steady state */
- external_flux =
- path->external_flux_term * source_get_power(green->scn->source, INF);
- }
+ power = green_path_power_contribution(green, path);
+ flux = green_path_flux_contribution(green, path);
+ external_flux = green_path_external_flux_contribution(green, path);
/* Compute path's end temperature */
switch(path->end_type) {
@@ -1279,6 +1374,33 @@ error:
}
res_T
+sdis_green_function_get_external_flux_terms_count
+ (const struct sdis_green_path* path_handle,
+ size_t* nterms)
+{
+ const struct green_path* path = NULL;
+ struct sdis_green_function* green = NULL;
+ res_T res = RES_OK;
+
+ if(!path_handle || !nterms) {
+ 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__;
+
+ *nterms = darray_extflux_terms_size_get(&path->extflux_terms);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
sdis_green_path_for_each_power_term
(struct sdis_green_path* path_handle,
sdis_process_medium_power_term_T func,
@@ -1354,25 +1476,41 @@ error:
}
res_T
-sdis_green_path_get_external_flux_term
+sdis_green_path_for_each_external_flux_terms
(struct sdis_green_path* path_handle,
- double* external_flux_term)
+ sdis_process_external_flux_terms_T func,
+ void* context)
{
const struct green_path* path = NULL;
struct sdis_green_function* green = NULL;
+ size_t i, n;
res_T res = RES_OK;
- if(!path_handle || !external_flux_term) {
+ if(!path_handle || !func) {
res = RES_BAD_ARG;
goto error;
}
- green = path_handle->green__; (void)green;
+ green = path_handle->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;
+ n = darray_extflux_terms_size_get(&path->extflux_terms);
+ if(n && !green->scn->source) {
+ /* In can't have external flux terms without an external source */
+ log_err(green->scn->dev, "%s: the external source is missing\n", FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ FOR_EACH(i, 0, n) {
+ const struct sdis_green_external_flux_terms* terms = NULL;
+ terms = darray_extflux_terms_cdata_get(&path->extflux_terms) + i;
+
+ res = func(green->scn->source, terms, context);
+ if(res != RES_OK) goto error;
+ }
exit:
return res;
@@ -1743,11 +1881,23 @@ error:
}
res_T
-green_path_add_external_flux_term
+green_path_add_external_flux_terms
(struct green_path_handle* handle,
- const double val) /* [W/m^2/sr] */
+ const struct sdis_green_external_flux_terms* terms)
{
- ASSERT(handle);
- handle->path->external_flux_term += val;
- return RES_OK;
+ res_T res = RES_OK;
+ ASSERT(handle && terms);
+
+ res = darray_extflux_terms_push_back(&handle->path->extflux_terms, terms);
+ if(res != RES_OK) {
+ log_err(handle->green->scn->dev,
+ "%s: cannot store external flux terms -- %s\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
}
diff --git a/src/sdis_green.h b/src/sdis_green.h
@@ -108,9 +108,9 @@ green_path_add_flux_term
const double term);
extern LOCAL_SYM res_T
-green_path_add_external_flux_term
+green_path_add_external_flux_terms
(struct green_path_handle* handle,
- const double term); /* [W/m^2/sr] */
+ const struct sdis_green_external_flux_terms* terms);
#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
@@ -50,6 +50,19 @@ struct brdf {
#define BRDF_NULL__ {0, 0}
static const struct brdf BRDF_NULL = BRDF_NULL__;
+/* Incident diffuse flux is made up of two components. One corresponds to the
+ * diffuse flux due to the reflection of the source on surfaces. The other is
+ * the diffuse flux due to the source's radiation scattering at least once in
+ * the environment. */
+struct incident_diffuse_flux {
+ double reflected; /* [W/m^2] */
+ double scattered; /* [W/m^2] */
+ double dir[3]; /* Direction along wich the scattered part was retrieved */
+};
+#define INCIDENT_DIFFUSE_FLUX_NULL__ {0, 0, {0,0,0}}
+static const struct incident_diffuse_flux INCIDENT_DIFFUSE_FLUX_NULL =
+ INCIDENT_DIFFUSE_FLUX_NULL__;
+
/* Reflect the V wrt the normal N. By convention V points outward the surface.
* In fact, this function is a double-precision version of the reflect_3d
* function. TODO Clean this "repeat" */
@@ -314,15 +327,14 @@ XD(compute_incident_diffuse_flux)
const double in_N[DIM], /* Surface normal. (Away from the surface) */
const double time,
const struct sXd(hit)* in_hit, /* Current intersection */
- double* out_flux) /* [W/m^2] */
+ struct incident_diffuse_flux* diffuse_flux) /* [W/m^2] */
{
struct sXd(hit) hit = SXD_HIT_NULL;
double pos[3] = {0}; /* In 3D for ray tracing ray to the source */
double dir[3] = {0}; /* Incident direction (toward the surface). Always 3D.*/
double N[3] = {0}; /* Surface normal. Always 3D */
- double incident_diffuse_flux = 0; /* [W/m^2] */
res_T res = RES_OK;
- ASSERT(in_pos && in_N && in_hit);
+ ASSERT(in_pos && in_N && in_hit && diffuse_flux);
/* Local copy of input argument */
dX(set)(pos, in_pos);
@@ -332,6 +344,8 @@ XD(compute_incident_diffuse_flux)
/* Sample a diffusive direction in 3D */
ssp_ran_hemisphere_cos(rng, N, dir, NULL);
+ *diffuse_flux = INCIDENT_DIFFUSE_FLUX_NULL;
+
for(;;) {
/* External sources */
struct source_sample src_sample = SOURCE_SAMPLE_NULL;
@@ -345,7 +359,7 @@ XD(compute_incident_diffuse_flux)
struct brdf_sample brdf_sample = BRDF_SAMPLE_NULL;
/* Miscellaneous */
- double L = 0; /* incident direct flux to bounce position */
+ double L = 0; /* incident flux to bounce position */
double wi[3] = {0}; /* Incident direction (outward the surface). Always 3D */
double vec[DIM] = {0}; /* Temporary variable */
@@ -353,7 +367,20 @@ XD(compute_incident_diffuse_flux)
/* Find the following surface along the direction of propagation */
XD(trace_ray)(scn, pos, dir, INF, &hit, &hit);
- if(SXD_HIT_NONE(&hit)) break; /* No surface */
+ if(SXD_HIT_NONE(&hit)) {
+ /* No surface. Handle the radiance emitted by the source and scattered at
+ * least once in the environment. Note that the value returned is not the
+ * actual scattered component of the incident diffuse flux: it relates
+ * to the radiance of the source scattered along the input dir at the
+ * given instant. It must therefore be multiplied by this radiance to
+ * obtain its real contribution. This trick makes it possible to manage
+ * the external flux in the green function. */
+ diffuse_flux->scattered = PI;
+ diffuse_flux->dir[0] = dir[0];
+ diffuse_flux->dir[1] = dir[1];
+ diffuse_flux->dir[2] = dir[2];
+ break;
+ }
/* Retrieve the current position and normal */
dX(add)(pos, pos, dX(muld)(vec, dir, hit.distance));
@@ -385,7 +412,7 @@ XD(compute_incident_diffuse_flux)
if(!SOURCE_SAMPLE_NONE(&src_sample)) {
const double Ld = XD(direct_contribution)(scn, &src_sample, pos, &hit);
- L = Ld; /* [W/m^2] */
+ L = Ld; /* [W/m^2/sr] */
}
/* Calculate the direct contribution of the rebound is diffuse */
@@ -401,21 +428,19 @@ XD(compute_incident_diffuse_flux)
/* The source is behind the surface */
if(cos_theta <= 0) {
- L = 0; /* [W/m^2] */
+ L = 0; /* [W/m^2/sr] */
/* The source is above the surface */
} else {
const double Ld = XD(direct_contribution)(scn, &src_sample, pos, &hit);
- L = Ld * cos_theta / (PI * src_sample.pdf); /* [W/m^2] */
+ L = Ld * cos_theta / (PI * src_sample.pdf); /* [W/m^2/sr] */
}
}
- incident_diffuse_flux += L;
+ diffuse_flux->reflected += L; /* [W/m^2/sr] */
}
-
- incident_diffuse_flux *= PI;
+ diffuse_flux->reflected *= PI; /* [W/m^2] */
exit:
- *out_flux = incident_diffuse_flux;
return res;
error:
goto exit;
@@ -431,15 +456,19 @@ XD(handle_external_net_flux)
const struct XD(handle_external_net_flux_args)* args,
struct XD(temperature)* T)
{
+ /* Terms to be registered in the green function */
+ struct sdis_green_external_flux_terms green =
+ SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL;
+
/* Sampling external sources */
struct source_sample src_sample = SOURCE_SAMPLE_NULL;
/* External flux */
+ struct incident_diffuse_flux incident_flux_diffuse = INCIDENT_DIFFUSE_FLUX_NULL;
double incident_flux = 0; /* [W/m^2] */
- 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] */
+ double net_flux_sc = 0; /* [W/m^2] */
/* Sampled path */
double N[3] = {0}; /* Normal. Always in 3D */
@@ -448,6 +477,7 @@ XD(handle_external_net_flux)
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
/* Miscellaneous */
+ double sum_h = 0;
double emissivity = 0; /* Emissivity */
double Ld = 0; /* Incident radiance [W/m^2/sr] */
double cos_theta = 0;
@@ -495,27 +525,50 @@ XD(handle_external_net_flux)
(scn, rng, frag.P, N, frag.time, args->hit, &incident_flux_diffuse);
if(res != RES_OK) goto error;
- /* Calculate the global incident flux */
- incident_flux = incident_flux_direct + incident_flux_diffuse; /* [W/m^2] */
-
- /* Calculate the net flux */
+ /* Calculate the incident flux without the part scattered by the environment.
+ * The latter depends on the source's diffuse radiance, not on its power. On
+ * the other hand, both the direct incident flux and the diffuse incident flux
+ * reflected by surfaces are linear with respect to the source power. This
+ * term can therefore be recorded in the green function in relation to this
+ * power, whereas the incident diffused flux coming from the scattered source
+ * radiance depends on the diffuse radiance of the source */
+ incident_flux = /* [W/m^2] */
+ 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] */
- /* 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 += external_flux_term * source_get_power(scn->source, frag.time);
+ /* Calculate the net flux from the radiance source scattered at least once by
+ * the medium */
+ net_flux_sc = incident_flux_diffuse.scattered * emissivity; /* [W/m^2] */
+
+ /* Until now, net flux has been calculated on the basis of source power and
+ * source diffuse radiance. What is actually calculated are the external flux
+ * terms of the green function. These must be multiplied by the source power
+ * and the source diffuse radiance, then added together to give the actual
+ * external flux */
+ sum_h = (args->h_radi + args->h_conv + args->h_cond);
+ green.term_wrt_power = net_flux / sum_h; /* [K/W] */
+ green.term_wrt_diffuse_radiance = net_flux_sc / sum_h; /* [K/W/m^2/sr] */
+ green.time = frag.time; /* [s] */
+ green.dir[0] = incident_flux_diffuse.dir[0];
+ green.dir[1] = incident_flux_diffuse.dir[1];
+ green.dir[2] = incident_flux_diffuse.dir[2];
+
+ T->value += green.term_wrt_power * source_get_power(scn->source, green.time);
+ if(green.term_wrt_diffuse_radiance) {
+ T->value +=
+ green.term_wrt_diffuse_radiance
+ * source_get_diffuse_radiance(scn->source, green.time, green.dir);
+ }
- /* Register the external net flux term */
+ /* Register the external net flux terms */
if(args->green_path) {
- res = green_path_add_external_flux_term(args->green_path, external_flux_term);
+ res = green_path_add_external_flux_terms(args->green_path, &green);
if(res != RES_OK) goto error;
}
diff --git a/src/sdis_source.c b/src/sdis_source.c
@@ -133,10 +133,18 @@ sdis_source_ref_put(struct sdis_source* src)
double
sdis_source_get_power(struct sdis_source* src, const double time /* [s] */)
{
- ASSERT(src);
return source_get_power(src, time);
}
+double
+sdis_source_get_diffuse_radiance
+ (struct sdis_source* src,
+ const double time, /* [s] */
+ const double dir[3])
+{
+ return source_get_diffuse_radiance(src, time, dir);
+}
+
unsigned
sdis_source_get_id(const struct sdis_source* source)
{
@@ -321,6 +329,20 @@ source_get_power(const struct sdis_source* src, const double time /* [s] */)
return src->spherical.power(time, src->spherical.data);
}
+double /* [W/perpendicular m^2/sr] */
+source_get_diffuse_radiance
+ (const struct sdis_source* src,
+ const double time /* [s] */,
+ const double dir[3])
+{
+ ASSERT(src);
+ if(src->spherical.diffuse_radiance == NULL) {
+ return 0;
+ } else {
+ return src->spherical.diffuse_radiance(time, dir, src->spherical.data);
+ }
+}
+
void
source_compute_signature(const struct sdis_source* src, hash256_T hash)
{
diff --git a/src/sdis_source_c.h b/src/sdis_source_c.h
@@ -66,6 +66,12 @@ source_get_power
(const struct sdis_source* source,
const double time); /* [s] */
+extern LOCAL_SYM double /* [W/m^2/sr] */
+source_get_diffuse_radiance
+ (const struct sdis_source* source,
+ const double time, /* [s] */
+ const double dir[3]);
+
extern LOCAL_SYM void
source_compute_signature
(const struct sdis_source* source,
diff --git a/src/test_sdis_draw_external_flux.c b/src/test_sdis_draw_external_flux.c
@@ -191,6 +191,17 @@ source_get_power(const double time, struct sdis_data* data)
return SOURCE_POWER; /* [W] */
}
+static double
+source_get_diffuse_radiance
+ (const double time,
+ const double dir[3],
+ struct sdis_data* data)
+{
+ (void)time, (void)data; /* Avoid the "unusued variable" warning */
+ CHK(d3_is_normalized(dir));
+ return 50;
+}
+
static struct sdis_source*
create_source(struct sdis_device* sdis)
{
@@ -199,6 +210,7 @@ create_source(struct sdis_device* sdis)
args.position = source_get_position;
args.power = source_get_power;
+ args.diffuse_radiance = source_get_diffuse_radiance;
args.data = NULL;
args.radius = 3e-1; /* [m] */
OK(sdis_spherical_source_create(sdis, &args, &source));
diff --git a/src/test_sdis_external_flux_with_diffuse_radiance.c b/src/test_sdis_external_flux_with_diffuse_radiance.c
@@ -0,0 +1,451 @@
+/* Copyright (C) 2016-2024 |Méso|Star> (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 "test_sdis_utils.h"
+#include "sdis.h"
+
+#include <rsys/rsys.h>
+#include <star/s3dut.h>
+
+/*
+ * The system is a sphere with an emissivity of 1. It is illuminated by an
+ * external spherical source that is sufficiently far away from the sphere
+ * radius to assume a constant flux density received per m^2 of surface
+ * perpendicular to the direction towards the source center.
+ *
+ * In such situation, the value of the surface temperature is equal to :
+ *
+ * Ts = [q/(4*sigma)+Tamb^4]^0.25
+ * q = P/(4PI*d^2)
+ *
+ * with Tamb the temperature of the environment, sigma the Boltzmann constant,
+ * P the source power and d the distance of the source from the center of
+ * the sphere.
+ *
+ * If one adds a diffuse radiance Ldiff, the surface temperature becomes :
+ *
+ * Ts = [(q/4+Ldiff*PI)/sigma+Tamb^4]^0.25
+ *
+ * The test checks that we retrieved these analatycal results by Monte Carlo
+ *
+ *
+ * R = 0.01 m External source
+ * __ d = 10 m ###
+ * T = ? -> / .\ ..............................##### r = 0.3 m
+ * \__/ #####
+ * E=1 ###
+ * P = 3^7 W
+ */
+
+/* The source */
+#define SOURCE_POWER 1.0e5 /* [W] */
+#define SOURCE_RADIUS 0.3 /* [m/fp_to_meter] */
+#define SOURCE_DISTANCE 10 /* [m/fp_to_meter] */
+
+/* Miscellaneous */
+#define SPHERE_RADIUS 0.01 /* [m/fp_to_meter] */
+#define SPHERE_T_REF 320 /* [K] */
+#define T_RAD 300.0 /* [K] */
+#define T_REF 300.0 /* [K] */
+#define T_RAD4 (T_RAD*T_RAD*T_RAD*T_RAD)
+
+#define NREALISATIONS 10000
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static struct s3dut_mesh*
+create_sphere(void)
+{
+ struct s3dut_mesh* mesh = NULL;
+ OK(s3dut_create_sphere(NULL, SPHERE_RADIUS, 128, 64, &mesh));
+ return mesh;
+}
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+#define MEDIUM_PROP(Type, Prop, Val) \
+ static double \
+ Type##_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+MEDIUM_PROP(medium, volumic_mass, 1) /* [kj/m^3] */
+MEDIUM_PROP(medium, calorific_capacity, 1) /* [J/K/Kg] */
+MEDIUM_PROP(medium, temperature, SDIS_TEMPERATURE_NONE) /* [K] */
+MEDIUM_PROP(solid, thermal_conductivity, 1) /* [W/m/K] */
+MEDIUM_PROP(solid, delta, (2*SPHERE_RADIUS)/10.0) /* [m] */
+#undef MEDIUM_PROP
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = medium_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = medium_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = medium_get_temperature;
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+static struct sdis_medium*
+create_fluid(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* fluid = NULL;
+
+ shader.calorific_capacity = medium_get_calorific_capacity;
+ shader.volumic_mass = medium_get_volumic_mass;
+ shader.temperature = medium_get_temperature;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &fluid));
+ return fluid;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
+{
+ (void)frag, (void)source_id, (void)data;
+ return 1;
+}
+
+static double
+interface_get_reference_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)frag, (void)data;
+ return SPHERE_T_REF; /* [K] */
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+
+ shader.front.emissivity = interface_get_emissivity;
+ shader.front.reference_temperature = interface_get_reference_temperature;
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * Radiative environment
+ ******************************************************************************/
+#define RADENV_PROP(Prop, Val) \
+ static double \
+ radenv_get_##Prop \
+ (const struct sdis_radiative_ray* ray, \
+ struct sdis_data* data) \
+ { \
+ (void)ray, (void)data; \
+ return (Val); /* [K] */ \
+ }
+RADENV_PROP(temperature, T_RAD)
+RADENV_PROP(reference_temperature, T_REF)
+#undef RADENV_PROP
+
+static struct sdis_radiative_env*
+create_radenv(struct sdis_device* sdis)
+{
+ struct sdis_radiative_env_shader shader = SDIS_RADIATIVE_ENV_SHADER_NULL;
+ struct sdis_radiative_env* radenv = NULL;
+
+ shader.temperature = radenv_get_temperature;
+ shader.reference_temperature = radenv_get_reference_temperature;
+ OK(sdis_radiative_env_create(sdis, &shader, NULL, &radenv));
+ return radenv;
+}
+
+/*******************************************************************************
+ * External source
+ ******************************************************************************/
+static void
+source_get_position
+ (const double time,
+ double pos[3],
+ struct sdis_data* data)
+{
+ (void)time, (void)data; /* Avoid the "unusued variable" warning */
+ pos[0] = SOURCE_DISTANCE;
+ pos[1] = 0;
+ pos[2] = 0;
+}
+
+static double
+source_get_power(const double time, struct sdis_data* data)
+{
+ (void)time, (void)data; /* Avoid the "unused variable" warning */
+ return SOURCE_POWER; /* [W] */
+}
+
+static double
+source_get_diffuse_radiance
+ (const double time,
+ const double dir[3],
+ struct sdis_data* data)
+{
+ const double* Ldiff = sdis_data_cget(data);;
+ (void)time, (void)dir; /* Avoid the "unused variable" warning */
+ return *Ldiff; /* [W/m^2/sr] */
+}
+
+static struct sdis_source*
+create_source(struct sdis_device* sdis, double** diffuse_radiance)
+{
+ struct sdis_spherical_source_create_args args =
+ SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_source* src = NULL;
+
+ OK(sdis_data_create(sdis, sizeof(double), ALIGNOF(double), NULL, &data));
+ *diffuse_radiance = sdis_data_get(data);
+ **diffuse_radiance = 0;
+
+ args.position = source_get_position;
+ args.power = source_get_power;
+ args.diffuse_radiance = source_get_diffuse_radiance;
+ args.data = data;
+ args.radius = SOURCE_RADIUS;
+ OK(sdis_spherical_source_create(sdis, &args, &src));
+ OK(sdis_data_ref_put(data));
+ return src;
+}
+
+/*******************************************************************************
+ * The scene
+ ******************************************************************************/
+struct scene_context {
+ const struct s3dut_mesh_data* mesh;
+ struct sdis_interface* interf;
+};
+static const struct scene_context SCENE_CONTEXT_NULL = {NULL, NULL};
+
+static void
+scene_get_indices(const size_t itri, size_t ids[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(ids && context && itri < context->mesh->nprimitives);
+ ids[0] = (unsigned)context->mesh->indices[itri*3+0];
+ ids[1] = (unsigned)context->mesh->indices[itri*3+1];
+ ids[2] = (unsigned)context->mesh->indices[itri*3+2];
+}
+
+static void
+scene_get_interface(const size_t itri, struct sdis_interface** interf, void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(interf && context && itri < context->mesh->nprimitives);
+ *interf = context->interf;
+}
+
+static void
+scene_get_position(const size_t ivert, double pos[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(pos && context && ivert < context->mesh->nvertices);
+ pos[0] = context->mesh->positions[ivert*3+0];
+ pos[1] = context->mesh->positions[ivert*3+1];
+ pos[2] = context->mesh->positions[ivert*3+2];
+}
+
+static struct sdis_scene*
+create_scene
+ (struct sdis_device* sdis,
+ const struct s3dut_mesh_data* mesh,
+ struct sdis_interface* interf,
+ struct sdis_source* source,
+ struct sdis_radiative_env* radenv)
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct scene_context context = SCENE_CONTEXT_NULL;
+
+ context.mesh = mesh;
+ context.interf = interf;
+
+ scn_args.get_indices = scene_get_indices;
+ scn_args.get_interface = scene_get_interface;
+ scn_args.get_position = scene_get_position;
+ scn_args.nprimitives = mesh->nprimitives;
+ scn_args.nvertices = mesh->nvertices;
+ scn_args.t_range[0] = MMIN(T_REF, SPHERE_T_REF);
+ scn_args.t_range[1] = MMAX(T_REF, SPHERE_T_REF);
+ scn_args.source = source;
+ scn_args.radenv = radenv;
+ scn_args.context = &context;
+ OK(sdis_scene_create(sdis, &scn_args, &scn));
+ return scn;
+}
+
+/*******************************************************************************
+ * Validations
+ ******************************************************************************/
+static double
+analytic_temperature(const double Ldiff)
+{
+ const double q = SOURCE_POWER/(4*PI*SOURCE_DISTANCE*SOURCE_DISTANCE);
+ const double Ts = pow((q/4 + Ldiff*PI)/BOLTZMANN_CONSTANT + T_RAD4, 0.25);
+ return Ts;
+}
+
+static void
+check
+ (struct sdis_scene* scn,
+ const struct s3dut_mesh_data* mesh,
+ const double Ldiff)
+{
+ struct sdis_solve_boundary_args args = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator = NULL;
+ enum sdis_side* sides = NULL;
+ double ref = 0;
+ size_t i = 0;
+
+ sides = mem_alloc(mesh->nprimitives*sizeof(*sides));
+ FOR_EACH(i, 0, mesh->nprimitives) sides[i] = SDIS_FRONT;
+
+ args.primitives = mesh->indices;
+ args.sides = sides;
+ args.nprimitives = mesh->nprimitives;
+ args.nrealisations = NREALISATIONS;
+ OK(sdis_solve_boundary(scn, &args, &estimator));
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ ref = analytic_temperature(Ldiff);
+ printf("Ts = %g ~ %g +/- %g\n", ref, T.E, T.SE);
+ OK(sdis_estimator_ref_put(estimator));
+
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+
+ mem_rm(sides);
+}
+
+static void
+solve_green(struct sdis_green_function* green, const double Ldiff)
+{
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator = NULL;
+ double ref = 0;
+
+ OK(sdis_green_function_solve(green, &estimator));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+
+ ref = analytic_temperature(Ldiff);
+ printf("Ts = %g ~ %g +/- %g\n", ref, T.E, T.SE);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+
+ OK(sdis_estimator_ref_put(estimator));
+}
+
+static void
+check_green
+ (struct sdis_scene* scn,
+ const struct s3dut_mesh_data* mesh,
+ double* Ldiff)
+{
+ struct sdis_solve_boundary_args args = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT;
+ struct sdis_green_function* green = NULL;
+ enum sdis_side* sides = NULL;
+ size_t i = 0;
+
+ CHK(Ldiff);
+
+ sides = mem_alloc(mesh->nprimitives*sizeof(*sides));
+ FOR_EACH(i, 0, mesh->nprimitives) sides[i] = SDIS_FRONT;
+
+ *Ldiff = 0;
+
+ args.primitives = mesh->indices;
+ args.sides = sides;
+ args.nprimitives = mesh->nprimitives;
+ args.nrealisations = NREALISATIONS;
+ OK(sdis_solve_boundary_green_function(scn, &args, &green));
+
+ solve_green(green, (*Ldiff = 50));
+ solve_green(green, (*Ldiff = 45));
+ solve_green(green, (*Ldiff = 40));
+
+ OK(sdis_green_function_ref_put(green));
+
+ mem_rm(sides);
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ /* Stardis */
+ struct sdis_device* dev = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_radiative_env* radenv = NULL;
+ struct sdis_scene* scene = NULL;
+ struct sdis_source* source = NULL;
+
+ /* Miscellaneous */
+ struct s3dut_mesh_data mesh;
+ struct s3dut_mesh* sphere = NULL;
+ double* Ldiff = NULL;
+
+ (void)argc, (void)argv;
+
+ OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev));
+
+ sphere = create_sphere();
+ OK(s3dut_mesh_get_data(sphere, &mesh));
+
+ fluid = create_fluid(dev);
+ solid = create_solid(dev);
+ interf = create_interface(dev, fluid, solid);
+ radenv = create_radenv(dev);
+ source = create_source(dev, &Ldiff);
+
+ scene = create_scene(dev, &mesh, interf, source, radenv);
+
+ check(scene, &mesh, (*Ldiff = 50));
+ check_green(scene, &mesh, Ldiff);
+
+ OK(s3dut_mesh_ref_put(sphere));
+ OK(sdis_device_ref_put(dev));
+ OK(sdis_interface_ref_put(interf));
+ OK(sdis_medium_ref_put(fluid));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_radiative_env_ref_put(radenv));
+ OK(sdis_scene_ref_put(scene));
+ OK(sdis_source_ref_put(source));
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -42,7 +42,7 @@ accum_power_terms(struct sdis_medium* mdm, const double power_term, void* ctx)
struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
struct sdis_data* data = NULL;
- double* power = ctx;
+ double* power = ctx; /* Power contribution [K] */
CHK(mdm && ctx);
CHK(sdis_medium_get_type(mdm) == SDIS_SOLID);
@@ -65,7 +65,7 @@ accum_flux_terms
struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
struct sdis_data* data = NULL;
- double* flux = ctx;
+ double* flux = ctx; /* Flux contribution [K] */
double phi;
CHK(interf && ctx);
@@ -84,6 +84,27 @@ accum_flux_terms
}
static res_T
+accum_extflux
+ (struct sdis_source* source,
+ const struct sdis_green_external_flux_terms* terms,
+ void* ctx)
+{
+ double* extflux = ctx; /* External flux contribution [K] */
+ double power = 0; /* [W] */
+ double diffuse_radiance = 0; /* [W/m^2/sr] */
+
+ CHK(source && terms && ctx);
+
+ power = sdis_source_get_power(source, terms->time);
+ diffuse_radiance = sdis_source_get_diffuse_radiance
+ (source, terms->time, terms->dir);
+
+ *extflux += terms->term_wrt_power * power;
+ *extflux += terms->term_wrt_diffuse_radiance * diffuse_radiance;
+ return RES_OK;
+}
+
+static res_T
solve_green_path(struct sdis_green_path* path, void* ctx)
{
struct sdis_green_path_end end = SDIS_GREEN_PATH_END_NULL;
@@ -99,13 +120,12 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
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;
double power = 0;
double flux = 0;
- double external_flux = 0; /* [W/m^2] */
+ double extflux = 0;
double time, temp = 0;
double weight = 0;
CHK(path && ctx);
@@ -130,23 +150,15 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
BA(sdis_green_path_for_each_flux_term(path, NULL, &acc));
OK(sdis_green_path_for_each_flux_term(path, accum_flux_terms, &flux));
+ BA(sdis_green_path_for_each_external_flux_terms(NULL, &accum_extflux, &extflux));
+ BA(sdis_green_path_for_each_external_flux_terms(path, NULL, &extflux));
+ OK(sdis_green_path_for_each_external_flux_terms(path, &accum_extflux, &extflux));
+
BA(sdis_green_path_get_elapsed_time(NULL, NULL));
BA(sdis_green_path_get_elapsed_time(path, NULL));
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 {
- /* NOTE: source power is assumed constant in time and is therefore retrieved
- * at steady state*/
- external_flux *= sdis_source_get_power(source, INF); /* [W] */
- }
-
BA(sdis_green_path_get_end(NULL, NULL));
BA(sdis_green_path_get_end(NULL, &end));
BA(sdis_green_path_get_end(path, NULL));
@@ -184,7 +196,7 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
default: FATAL("Unreachable code.\n"); break;
}
- weight = temp + power + external_flux + flux;
+ weight = temp + power + extflux + flux;
acc->sum += weight;
acc->sum2 += weight*weight;
