commit 35309912f04233443f9dbe897286cb22fc9c4df2
parent def3a4cdde0f1e5f99b567edd8ba7643a960afdd
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 15 Apr 2024 08:35:22 +0200
Merge branch 'feature_per_source_rad_props' into develop
Diffstat:
41 files changed, 2124 insertions(+), 1331 deletions(-)
diff --git a/Makefile b/Makefile
@@ -43,6 +43,7 @@ SRC =\
src/sdis_log.c\
src/sdis_medium.c\
src/sdis_misc.c\
+ src/sdis_radiative_env.c\
src/sdis_realisation.c\
src/sdis_scene.c\
src/sdis_solve.c\
@@ -188,6 +189,7 @@ TEST_SRC =\
src/test_sdis_interface.c\
src/test_sdis_medium.c\
src/test_sdis_picard.c\
+ src/test_sdis_radiative_env.c\
src/test_sdis_scene.c\
src/test_sdis_solid_random_walk_robustness.c\
src/test_sdis_solve_probe.c\
@@ -197,11 +199,11 @@ TEST_SRC =\
src/test_sdis_solve_probe3_2d.c\
src/test_sdis_source.c\
src/test_sdis_transcient.c\
- src/test_sdis_unstationary_atm.c\
src/test_sdis_unsteady.c\
src/test_sdis_unsteady_1d.c\
src/test_sdis_unsteady_analytic_profile.c\
src/test_sdis_unsteady_analytic_profile_2d.c\
+ src/test_sdis_unsteady_atm.c\
src/test_sdis_volumic_power.c\
src/test_sdis_volumic_power4.c
TEST_SRC_LONG =\
@@ -297,16 +299,17 @@ src/test_sdis_flux_with_h.d \
src/test_sdis_interface.d \
src/test_sdis_medium.d \
src/test_sdis_picard.d \
+src/test_sdis_radiative_env.d \
src/test_sdis_solve_probe.d \
src/test_sdis_solve_probe_2d.d \
src/test_sdis_solve_probe2_2d.d \
src/test_sdis_solve_probe3_2d \
src/test_sdis_source.d \
src/test_sdis_transcient.d \
-src/test_sdis_unstationary_atm.d \
src/test_sdis_unsteady.d \
src/test_sdis_unsteady_1d.d \
src/test_sdis_unsteady_analytic_profile_2d.d \
+src/test_sdis_unsteady_atm.d \
src/test_sdis_utils.d \
src/test_sdis_volumic_power.d \
src/test_sdis_volumic_power2.d \
@@ -331,16 +334,17 @@ src/test_sdis_flux_with_h.o \
src/test_sdis_interface.o \
src/test_sdis_medium.o \
src/test_sdis_picard.o \
+src/test_sdis_radiative_env.o \
src/test_sdis_solve_probe.o \
src/test_sdis_solve_probe_2d.o \
src/test_sdis_solve_probe2_2d.o \
src/test_sdis_solve_probe3_2d.o \
src/test_sdis_source.o \
src/test_sdis_transcient.o \
-src/test_sdis_unstationary_atm.o \
src/test_sdis_unsteady.o \
src/test_sdis_unsteady_1d.o \
src/test_sdis_unsteady_analytic_profile_2d.o \
+src/test_sdis_unsteady_atm.o \
src/test_sdis_utils.o \
src/test_sdis_volumic_power.o \
src/test_sdis_volumic_power2.o \
@@ -365,16 +369,17 @@ test_sdis_flux_with_h \
test_sdis_interface \
test_sdis_medium \
test_sdis_picard \
+test_sdis_radiative_env \
test_sdis_solve_probe \
test_sdis_solve_probe_2d \
test_sdis_solve_probe2_2d \
test_sdis_solve_probe3_2d \
test_sdis_source \
test_sdis_transcient \
-test_sdis_unstationary_atm \
test_sdis_unsteady \
test_sdis_unsteady_1d \
test_sdis_unsteady_analytic_profile_2d \
+test_sdis_unsteady_atm \
test_sdis_volumic_power \
test_sdis_volumic_power2 \
test_sdis_volumic_power2_2d \
diff --git a/src/sdis.h b/src/sdis.h
@@ -22,7 +22,9 @@
#include <rsys/hash.h>
#include <rsys/rsys.h>
-#include <float.h>
+
+#include <float.h> /* DBL_MAX */
+#include <limits.h> /* UINT_MAX */
/* Library symbol management */
#if defined(SDIS_SHARED_BUILD)
@@ -55,6 +57,9 @@
#define SDIS_TEMPERATURE_IS_KNOWN(Temp) (!IS_NaN(Temp))
#define SDIS_TEMPERATURE_IS_UNKNOWN(Temp) (IS_NaN(Temp))
+/* Identifier of the internal source of radiation */
+#define SDIS_INTERN_SOURCE_ID UINT_MAX
+
/* Forward declaration of external opaque data types */
struct logger;
struct mem_allocator;
@@ -75,6 +80,7 @@ struct sdis_estimator_buffer;
struct sdis_green_function;
struct sdis_interface;
struct sdis_medium;
+struct sdis_radiative_env; /* Radiative environment */
struct sdis_scene;
struct sdis_source;
@@ -103,7 +109,7 @@ enum sdis_diffusion_algorithm {
SDIS_DIFFUSION_NONE = SDIS_DIFFUSION_ALGORITHMS_COUNT__
};
-/* Random walk vertex, i.e. a spatiotemporal position at a given step of the
+/* Random walk vertex, i.e. a spatio-temporal position at a given step of the
* random walk. */
struct sdis_rwalk_vertex {
double P[3]; /* World space position */
@@ -113,7 +119,7 @@ struct sdis_rwalk_vertex {
static const struct sdis_rwalk_vertex SDIS_RWALK_VERTEX_NULL =
SDIS_RWALK_VERTEX_NULL__;
-/* Spatiotemporal position onto an interface. As a random walk vertex, it
+/* Spatio-temporal position onto an interface. As a random walk vertex, it
* stores the position and time of the random walk, but since it lies onto an
* interface, it has additionnal parameters as the normal of the interface and
* the parametric coordinate of the position onto the interface */
@@ -128,6 +134,14 @@ struct sdis_interface_fragment {
static const struct sdis_interface_fragment SDIS_INTERFACE_FRAGMENT_NULL =
SDIS_INTERFACE_FRAGMENT_NULL__;
+/* Ray traced in radiative environment */
+struct sdis_radiative_ray {
+ double dir[3]; /* Direction */
+};
+#define SDIS_RADIATIVE_RAY_NULL__ {{0,0,0}}
+static const struct sdis_radiative_ray SDIS_RADIATIVE_RAY_NULL=
+ SDIS_RADIATIVE_RAY_NULL__;
+
/* Input arguments of the sdis_device_create function */
struct sdis_device_create_args {
struct logger* logger; /* NULL <=> default logger */
@@ -234,6 +248,20 @@ typedef double
(const struct sdis_interface_fragment* frag, /* Interface position */
struct sdis_data* data); /* User data */
+/* Type of functor for obtaining the spatio temporal physical properties of an
+ * interface, as a function of the radiation source */
+typedef double
+(*sdis_radiative_interface_getter_T)
+ (const struct sdis_interface_fragment* frag, /* Interface position */
+ const unsigned source_id, /* Identifier of the radiation source */
+ struct sdis_data* data); /* User data */
+
+/* Type of functor for obtaining radiative environment properties */
+typedef double
+(*sdis_radiative_ray_getter_T)
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data);
+
/* Define the physical properties of a solid */
struct sdis_solid_shader {
/* Properties */
@@ -290,8 +318,8 @@ struct sdis_interface_side_shader {
/* Control the emissivity of the interface. May be NULL for solid/solid
* interface or if the emissivity is 0 onto the whole interface. */
- sdis_interface_getter_T emissivity; /* Overall emissivity. */
- sdis_interface_getter_T specular_fraction; /* Specular part in [0,1] */
+ sdis_radiative_interface_getter_T emissivity; /* Overall emissivity */
+ sdis_radiative_interface_getter_T specular_fraction; /* Specular part in [0,1] */
/* Reference temperature used in Picard 1 */
sdis_interface_getter_T reference_temperature;
@@ -326,6 +354,14 @@ struct sdis_interface_shader {
static const struct sdis_interface_shader SDIS_INTERFACE_SHADER_NULL =
SDIS_INTERFACE_SHADER_NULL__;
+struct sdis_radiative_env_shader {
+ sdis_radiative_ray_getter_T temperature; /* [K] */
+ sdis_radiative_ray_getter_T reference_temperature; /* [K] */
+};
+#define SDIS_RADIATIVE_ENV_SHADER_NULL__ {NULL, NULL}
+static const struct sdis_radiative_env_shader SDIS_RADIATIVE_ENV_SHADER_NULL =
+ SDIS_RADIATIVE_ENV_SHADER_NULL__;
+
/*******************************************************************************
* Registered heat path data types
******************************************************************************/
@@ -371,35 +407,39 @@ typedef res_T
******************************************************************************/
enum sdis_green_path_end_type {
SDIS_GREEN_PATH_END_AT_INTERFACE,
+ SDIS_GREEN_PATH_END_AT_RADIATIVE_ENV,
SDIS_GREEN_PATH_END_IN_VOLUME,
- SDIS_GREEN_PATH_END_RADIATIVE,
SDIS_GREEN_PATH_END_TYPES_COUNT__,
SDIS_GREEN_PATH_END_ERROR = SDIS_GREEN_PATH_END_TYPES_COUNT__
};
-enum sdis_point_type {
- SDIS_FRAGMENT,
- SDIS_VERTEX,
- SDIS_POINT_TYPES_COUNT__,
- SDIS_POINT_NONE = SDIS_POINT_TYPES_COUNT__
-};
-
/* Spatio temporal point */
-struct sdis_point {
+struct sdis_green_path_end {
union {
+ /* Path end in volume */
struct {
struct sdis_medium* medium;
struct sdis_rwalk_vertex vertex;
- } mdmvert; /* Medium and a vertex into it */
+ } mdmvert;
+ /* Path end at interface */
struct {
struct sdis_interface* intface;
struct sdis_interface_fragment fragment;
- } itfrag; /* Interface and a fragmetn onto it */
+ } itfrag;
+ /* Path end in radiative environement */
+ struct {
+ struct sdis_radiative_env* radenv;
+ struct sdis_radiative_ray ray;
+ } radenvray;
} data;
- enum sdis_point_type type;
+ enum sdis_green_path_end_type type;
};
-#define SDIS_POINT_NULL__ { {{NULL, SDIS_RWALK_VERTEX_NULL__}}, SDIS_POINT_NONE}
-static const struct sdis_point SDIS_POINT_NULL = SDIS_POINT_NULL__;
+#define SDIS_GREEN_PATH_END_NULL__ { \
+ {{NULL, SDIS_RWALK_VERTEX_NULL__}}, \
+ SDIS_GREEN_PATH_END_ERROR \
+}
+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 */
typedef res_T
@@ -450,18 +490,6 @@ typedef void
double pos[], /* Output list of vertex coordinates */
void* ctx);
-struct sdis_ambient_radiative_temperature {
- double temperature; /* In Kelvin */
- double reference; /* Used to linearise the radiative transfer */
-};
-#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__;
-
struct sdis_scene_create_args {
/* Functors to retrieve the geometric description */
sdis_get_primitive_indices_T get_indices;
@@ -474,7 +502,6 @@ struct sdis_scene_create_args {
size_t nprimitives; /* #primitives, i.e. #segments or #triangles */
size_t nvertices; /* #vertices */
double fp_to_meter; /* Scale factor used to convert a float in meter */
- struct sdis_ambient_radiative_temperature trad; /* Ambient radiative temp */
/* Min/max temperature used to linearise the radiative temperature */
double t_range[2];
@@ -482,6 +509,10 @@ struct sdis_scene_create_args {
/* External source. Can be NULL <=> no external flux will be calculated on
* scene interfaces */
struct sdis_source* source;
+
+ /* Radiative environment. Can be NULL <=> sampled radiative trajectories
+ * cannot (in fact must not) reach the surrounding environment */
+ struct sdis_radiative_env* radenv;
};
#define SDIS_SCENE_CREATE_ARGS_DEFAULT__ { \
@@ -492,9 +523,9 @@ struct sdis_scene_create_args {
0, /* #primitives */ \
0, /* #vertices */ \
1.0, /* #Floating point to meter scale factor */ \
- SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__,/* Ambient radiative temperature */\
{SDIS_TEMPERATURE_NONE, SDIS_TEMPERATURE_NONE}, /* Temperature range */ \
- NULL /* source */ \
+ NULL, /* source */ \
+ NULL /* Radiative environement */ \
}
static const struct sdis_scene_create_args SDIS_SCENE_CREATE_ARGS_DEFAULT =
SDIS_SCENE_CREATE_ARGS_DEFAULT__;
@@ -1042,6 +1073,34 @@ sdis_interface_get_id
(const struct sdis_interface* interf);
/*******************************************************************************
+ * API of the radiative environment. Describes the system when the sampled
+ * radiative paths reach infinity.
+ ******************************************************************************/
+SDIS_API res_T
+sdis_radiative_env_create
+ (struct sdis_device* dev,
+ const struct sdis_radiative_env_shader* shader,
+ struct sdis_data* data, /* Data sent to the shader. May be NULL */
+ struct sdis_radiative_env** radenv);
+
+SDIS_API res_T
+sdis_radiative_env_ref_get
+ (struct sdis_radiative_env* radenv);
+
+SDIS_API res_T
+sdis_radiative_env_ref_put
+ (struct sdis_radiative_env* radenv);
+
+SDIS_API res_T
+sdis_radiative_env_get_shader
+ (struct sdis_radiative_env* radenv,
+ struct sdis_radiative_env_shader* shader);
+
+SDIS_API struct sdis_data*
+sdis_radiative_env_get_data
+ (struct sdis_radiative_env* radenv);
+
+/*******************************************************************************
* External source API. When a scene has external sources, an external flux
* (in both its direct and diffuse parts) is imposed on the interfaces.
******************************************************************************/
@@ -1064,6 +1123,10 @@ sdis_source_get_power
(struct sdis_source* source,
const double time); /* [s] */
+SDIS_API unsigned
+sdis_source_get_id
+ (const struct sdis_source* source);
+
/*******************************************************************************
* A scene is a collection of primitives. Each primitive is the geometric
* support of the interface between 2 media.
@@ -1137,19 +1200,6 @@ sdis_scene_set_fp_to_meter
(struct sdis_scene* scn,
const double fp_to_meter);
-/* Get scene's ambient radiative temperature */
-SDIS_API res_T
-sdis_scene_get_ambient_radiative_temperature
- (const struct sdis_scene* scn,
- struct sdis_ambient_radiative_temperature* trad);
-
-/* Set scene's ambient radiative temperature. If set negative, any sample
- * ending in ambient radiative temperature will fail */
-SDIS_API res_T
-sdis_scene_set_ambient_radiative_temperature
- (struct sdis_scene* scn,
- const struct sdis_ambient_radiative_temperature* trad);
-
/* Get scene's minimum/maximum temperature */
SDIS_API res_T
sdis_scene_get_temperature_range
@@ -1255,6 +1305,12 @@ sdis_scene_get_source
(struct sdis_scene* scn,
struct sdis_source** src); /* The returned pointer can be NULL <=> no source */
+SDIS_API res_T
+sdis_scene_get_radiative_env
+ (struct sdis_scene* scn,
+ /* The returned pointer can be NULL, i.e. there is no radiative environement*/
+ struct sdis_radiative_env** radenv);
+
/*******************************************************************************
* An estimator stores the state of a simulation
******************************************************************************/
@@ -1405,18 +1461,12 @@ sdis_green_path_get_elapsed_time
(struct sdis_green_path* path_handle,
double* elapsed);
-/* Retrieve the path's end type. */
-SDIS_API res_T
-sdis_green_path_get_end_type
- (struct sdis_green_path* path,
- enum sdis_green_path_end_type* type);
-
-/* Retrieve the spatio-temporal end point of a path used to estimate the green
- * function. Return RES_BAD_OP for paths ending radiative. */
+/* Retrieve the spatio-temporal limit of a path used to estimate the green
+ * function */
SDIS_API res_T
-sdis_green_path_get_limit_point
+sdis_green_path_get_end
(struct sdis_green_path* path,
- struct sdis_point* pt);
+ struct sdis_green_path_end* end);
/* Retrieve the green function the path belongs to */
SDIS_API res_T
diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h
@@ -93,11 +93,15 @@ struct XD(rwalk) {
struct sdis_rwalk_vertex vtx; /* Position and time of the Random walk */
struct sdis_medium* mdm; /* Medium in which the random walk lies */
struct sXd(hit) hit; /* Hit of the random walk */
+
+ /* Direction along which the random walk reached the radiative environment */
+ double dir[3];
+
double elapsed_time;
enum sdis_side hit_side;
};
static const struct XD(rwalk) XD(RWALK_NULL) = {
- SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__, 0, SDIS_SIDE_NULL__
+ SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__, {0,0,0}, 0, SDIS_SIDE_NULL__
};
struct XD(temperature) {
diff --git a/src/sdis_device.c b/src/sdis_device.c
@@ -316,8 +316,10 @@ device_release(ref_T* ref)
if(dev->logger == &dev->logger__) logger_release(&dev->logger__);
ASSERT(flist_name_is_empty(&dev->interfaces_names));
ASSERT(flist_name_is_empty(&dev->media_names));
+ ASSERT(flist_name_is_empty(&dev->source_names));
flist_name_release(&dev->interfaces_names);
flist_name_release(&dev->media_names);
+ flist_name_release(&dev->source_names);
#ifdef SDIS_ENABLE_MPI
if(dev->mpi_mutex) mutex_destroy(dev->mpi_mutex);
str_release(&dev->mpi_err_str);
@@ -366,6 +368,7 @@ sdis_device_create
ref_init(&dev->ref);
flist_name_init(allocator, &dev->interfaces_names);
flist_name_init(allocator, &dev->media_names);
+ flist_name_init(allocator, &dev->source_names);
#ifdef SDIS_ENABLE_MPI
str_init(allocator, &dev->mpi_err_str);
#endif
diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h
@@ -61,6 +61,7 @@ struct sdis_device {
struct flist_name interfaces_names;
struct flist_name media_names;
+ struct flist_name source_names;
struct s2d_device* s2d_dev;
struct s3d_device* s3d_dev;
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -20,6 +20,7 @@
#include "sdis_log.h"
#include "sdis_medium_c.h"
#include "sdis_misc.h"
+#include "sdis_radiative_env_c.h"
#include "sdis_scene_c.h"
#include "sdis_source_c.h"
@@ -90,6 +91,7 @@ struct green_path {
union {
struct sdis_rwalk_vertex vertex;
struct sdis_interface_fragment fragment;
+ struct sdis_radiative_ray ray;
} limit;
unsigned limit_id; /* Identifier of the limit medium/interface */
enum sdis_green_path_end_type end_type;
@@ -110,6 +112,7 @@ green_path_init(struct mem_allocator* allocator, struct green_path* path)
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;
path->limit_id = UINT_MAX;
path->end_type = SDIS_GREEN_PATH_END_TYPES_COUNT__;
path->ilast_medium = UINT16_MAX;
@@ -405,8 +408,6 @@ green_function_solve_path
const size_t ipath,
double* weight)
{
- 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;
const struct green_path* path = NULL;
@@ -415,6 +416,7 @@ green_function_solve_path
struct sdis_scene* scn = NULL;
struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+ struct sdis_radiative_ray ray = SDIS_RADIATIVE_RAY_NULL;
double power;
double flux;
double external_flux;
@@ -466,23 +468,26 @@ green_function_solve_path
frag = path->limit.fragment;
end_temperature = interface_side_get_temperature(interf, &frag);
break;
+ case SDIS_GREEN_PATH_END_AT_RADIATIVE_ENV:
+ SDIS(green_function_get_scene(green, &scn));
+ ray = path->limit.ray;
+ end_temperature = radiative_env_get_temperature(green->scn->radenv, &ray);
+ break;
case SDIS_GREEN_PATH_END_IN_VOLUME:
medium = green_function_fetch_medium(green, path->limit_id);
vtx = path->limit.vertex;
end_temperature = medium_get_temperature(medium, &vtx);
break;
- case SDIS_GREEN_PATH_END_RADIATIVE:
- SDIS(green_function_get_scene(green, &scn));
- SDIS(scene_get_ambient_radiative_temperature(scn, &trad));
- end_temperature = trad.temperature;
- if(end_temperature < 0) { /* Cannot be negative if used */
- res = RES_BAD_ARG;
- goto error;
- }
- break;
default: FATAL("Unreachable code.\n"); break;
}
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(end_temperature)) {
+ log_err(green->scn->dev,
+ "%s: unknown boundary/initial condition.\n", FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
/* Compute the path weight */
*weight = power + flux + external_flux + end_temperature;
@@ -1149,39 +1154,15 @@ error:
}
res_T
-sdis_green_path_get_end_type
- (struct sdis_green_path* path_handle, enum sdis_green_path_end_type* type)
-{
- const struct green_path* path = NULL;
- struct sdis_green_function* green = NULL;
- res_T res = RES_OK;
-
- if(!path_handle || !type) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- 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__;
- *type = path->end_type;
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-res_T
-sdis_green_path_get_limit_point
- (struct sdis_green_path* path_handle, struct sdis_point* pt)
+sdis_green_path_get_end
+ (struct sdis_green_path* path_handle,
+ struct sdis_green_path_end* end)
{
const struct green_path* path = NULL;
struct sdis_green_function* green = NULL;
res_T res = RES_OK;
- if(!path_handle || !pt) {
+ if(!path_handle || !end) {
res = RES_BAD_ARG;
goto error;
}
@@ -1190,19 +1171,21 @@ sdis_green_path_get_limit_point
ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths));
path = darray_green_path_cdata_get(&green->paths) + path_handle->id__;
+ end->type = path->end_type;
switch(path->end_type) {
case SDIS_GREEN_PATH_END_AT_INTERFACE:
- pt->data.itfrag.intface = green_function_fetch_interf(green, path->limit_id);
- pt->data.itfrag.fragment = path->limit.fragment;
- pt->type = SDIS_FRAGMENT;
+ end->data.itfrag.intface = green_function_fetch_interf(green, path->limit_id);
+ end->data.itfrag.fragment = path->limit.fragment;
+ break;
+ case SDIS_GREEN_PATH_END_AT_RADIATIVE_ENV:
+ end->data.radenvray.radenv = green->scn->radenv;
+ end->data.radenvray.ray = path->limit.ray;
break;
case SDIS_GREEN_PATH_END_IN_VOLUME:
- pt->data.mdmvert.medium = green_function_fetch_medium(green, path->limit_id);
- pt->data.mdmvert.vertex = path->limit.vertex;
- pt->type = SDIS_VERTEX;
+ end->data.mdmvert.medium = green_function_fetch_medium(green, path->limit_id);
+ end->data.mdmvert.vertex = path->limit.vertex;
break;
- case SDIS_GREEN_PATH_END_RADIATIVE:
case SDIS_GREEN_PATH_END_ERROR:
res = RES_BAD_OP;
goto error;
@@ -1620,14 +1603,16 @@ green_path_set_limit_vertex
}
res_T
-green_path_set_limit_radiative
+green_path_set_limit_radiative_ray
(struct green_path_handle* handle,
+ const struct sdis_radiative_ray* ray,
const double elapsed_time)
{
ASSERT(handle);
ASSERT(handle->path->end_type == SDIS_GREEN_PATH_END_TYPES_COUNT__);
handle->path->elapsed_time = elapsed_time;
- handle->path->end_type = SDIS_GREEN_PATH_END_RADIATIVE;
+ handle->path->limit.ray = *ray;
+ handle->path->end_type = SDIS_GREEN_PATH_END_AT_RADIATIVE_ENV;
return RES_OK;
}
diff --git a/src/sdis_green.h b/src/sdis_green.h
@@ -84,8 +84,9 @@ green_path_set_limit_vertex
const double elapsed_time);
extern LOCAL_SYM res_T
-green_path_set_limit_radiative
+green_path_set_limit_radiative_ray
(struct green_path_handle* handle,
+ const struct sdis_radiative_ray* ray,
const double elapsed_time);
extern LOCAL_SYM res_T
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
@@ -272,23 +272,27 @@ XD(setup_fragment)
static INLINE res_T
XD(setup_brdf)
(struct sdis_device* dev,
+ const struct sdis_source* src,
struct brdf* brdf,
const struct sdis_interface* interf,
const struct sdis_interface_fragment* frag)
{
double epsilon = 0;
double alpha = 0;
+ unsigned src_id = 0;
res_T res = RES_OK;
ASSERT(brdf && frag);
ASSERT((frag->side == SDIS_FRONT
&& sdis_medium_get_type(interf->medium_front) == SDIS_FLUID)
|| sdis_medium_get_type(interf->medium_back) == SDIS_FLUID);
- epsilon = interface_side_get_emissivity(interf, frag);
+ src_id = sdis_source_get_id(src);
+
+ epsilon = interface_side_get_emissivity(interf, src_id, frag);
res = interface_side_check_emissivity(dev, epsilon, frag->P, frag->time);
if(res != RES_OK) goto error;
- alpha = interface_side_get_specular_fraction(interf, frag);
+ alpha = interface_side_get_specular_fraction(interf, src_id, frag);
res = interface_side_check_specular_fraction(dev, alpha, frag->P, frag->time);
if(res != RES_OK) goto error;
@@ -364,7 +368,7 @@ XD(compute_incident_diffuse_flux)
res = check_interface(interf, &frag);
if(res != RES_OK) goto error;
- XD(setup_brdf)(scn->dev, &brdf, interf, &frag);
+ XD(setup_brdf)(scn->dev, scn->source, &brdf, interf, &frag);
/* Check if path is absorbed */
if(ssp_rng_canonical(rng) < brdf.emissivity) break;
@@ -447,6 +451,7 @@ XD(handle_external_net_flux)
double emissivity = 0; /* Emissivity */
double Ld = 0; /* Incident radiance [W/m^2/sr] */
double cos_theta = 0;
+ unsigned src_id = 0;
int handle_flux = 0;
res_T res = RES_OK;
ASSERT(scn && args && T);
@@ -494,7 +499,8 @@ XD(handle_external_net_flux)
incident_flux = incident_flux_direct + incident_flux_diffuse; /* [W/m^2] */
/* Calculate the net flux */
- emissivity = interface_side_get_emissivity(args->interf, &frag);
+ 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] */
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h
@@ -18,6 +18,7 @@
#include "sdis_interface_c.h"
#include "sdis_medium_c.h"
#include "sdis_misc.h"
+#include "sdis_radiative_env_c.h"
#include "sdis_scene_c.h"
#include <star/ssp.h>
@@ -43,7 +44,11 @@ XD(rwalk_get_Tref)
* fetches the ambient radiative temperature. We do not use the limit
* conditions as the reference temperature to make the sampled paths
* independant of them. */
- Tref = scn->trad.reference;
+ struct sdis_radiative_ray ray = SDIS_RADIATIVE_RAY_NULL;
+ ray.dir[0] = rwalk->dir[0];
+ ray.dir[1] = rwalk->dir[1];
+ ray.dir[2] = rwalk->dir[2];
+ Tref = radiative_env_get_reference_temperature(scn->radenv, &ray);
} else {
struct sdis_interface_fragment frag;
struct sdis_interface* interf = NULL;
@@ -154,7 +159,8 @@ XD(solid_fluid_boundary_picard1_path)
delta = solid_get_delta(solid, &rwalk->vtx);
/* Fetch the boundary emissivity */
- epsilon = interface_side_get_emissivity(interf, &frag_fluid);
+ epsilon = interface_side_get_emissivity
+ (interf, SDIS_INTERN_SOURCE_ID, &frag_fluid);
if(epsilon <= 0) {
Tref = 0;
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h
@@ -209,7 +209,8 @@ XD(solid_fluid_boundary_picardN_path)
delta = solid_get_delta(solid, &rwalk->vtx);
/* Fetch the boundary emissivity */
- epsilon = interface_side_get_emissivity(interf, &frag_fluid);
+ epsilon = interface_side_get_emissivity
+ (interf, SDIS_INTERN_SOURCE_ID, &frag_fluid);
/* Note that the reinjection distance is *FIXED*. It MUST ensure that the
* orthogonal distance from the boundary to the reinjection point is at most
diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h
@@ -20,6 +20,7 @@
#include "sdis_log.h"
#include "sdis_medium_c.h"
#include "sdis_misc.h"
+#include "sdis_radiative_env_c.h"
#include "sdis_scene_c.h"
#include <star/ssp.h>
@@ -77,44 +78,52 @@ XD(trace_radiative_path)
(scn->sXd(view), pos, dir, range, &filter_data, &rwalk->hit));
#endif
if(SXD_HIT_NONE(&rwalk->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__;
- if(SDIS_TEMPERATURE_IS_KNOWN(scn->trad.temperature)) {
- T->value += scn->trad.temperature;
- T->done = 1;
-
- if(ctx->green_path) {
- res = green_path_set_limit_radiative
- (ctx->green_path, rwalk->elapsed_time);
- if(res != RES_OK) goto error;
- }
- if(ctx->heat_path) {
- const float empirical_dst = 0.1f;
- struct sdis_rwalk_vertex vtx;
-
-
- 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;
- }
- break;
- } else {
+ 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 reaches an invalid ambient radiative temperature "
- "of `%gK' at position `%g %g %g'. This may be due to numerical "
- "inaccuracies or to inconsistency in the simulated system (eg: "
- "unclosed geometry). For systems where the random walks can reach "
- "such temperature, one has to setup a valid ambient radiative "
- "temperature, i.e. it must be greater or equal to 0.\n",
- FUNC_NAME,
- scn->trad.temperature,
- SPLIT3(rwalk->vtx.P));
+ "%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;
}
/* Define the hit side */
@@ -134,7 +143,7 @@ XD(trace_radiative_path)
XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
/* Fetch the interface emissivity */
- epsilon = interface_side_get_emissivity(interf, &frag);
+ epsilon = interface_side_get_emissivity(interf, SDIS_INTERN_SOURCE_ID, &frag);
if(epsilon > 1 || epsilon < 0) {
log_err(scn->dev,
"%s: invalid overall emissivity `%g' at position `%g %g %g'.\n",
@@ -176,7 +185,7 @@ XD(trace_radiative_path)
goto error;
}
}
- alpha = interface_side_get_specular_fraction(interf, &frag);
+ alpha = interface_side_get_specular_fraction(interf, SDIS_INTERN_SOURCE_ID, &frag);
r = ssp_rng_canonical(rng);
if(r < alpha) { /* Sample specular part */
reflect_3d(dir, f3_minus(dir, dir), N);
diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h
@@ -143,6 +143,7 @@ interface_side_get_flux
static INLINE double
interface_side_get_emissivity
(const struct sdis_interface* interf,
+ const unsigned source_id,
const struct sdis_interface_fragment* frag)
{
const struct sdis_interface_side_shader* shader;
@@ -152,12 +153,15 @@ interface_side_get_emissivity
case SDIS_BACK: shader = &interf->shader.back; break;
default: FATAL("Unreachable code\n"); break;
}
- return shader->emissivity ? shader->emissivity(frag, interf->data) : 0;
+ return shader->emissivity
+ ? shader->emissivity(frag, source_id, interf->data)
+ : 0;
}
static INLINE double
interface_side_get_specular_fraction
(const struct sdis_interface* interf,
+ const unsigned source_id,
const struct sdis_interface_fragment* frag)
{
const struct sdis_interface_side_shader* shader;
@@ -168,7 +172,8 @@ interface_side_get_specular_fraction
default: FATAL("Unreachable code\n"); break;
}
return shader->specular_fraction
- ? shader->specular_fraction(frag, interf->data) : 0;
+ ? shader->specular_fraction(frag, source_id, interf->data)
+ : 0;
}
static INLINE double
diff --git a/src/sdis_radiative_env.c b/src/sdis_radiative_env.c
@@ -0,0 +1,108 @@
+/* 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 "sdis_device_c.h"
+#include "sdis_log.h"
+#include "sdis_radiative_env_c.h"
+
+#include <rsys/mem_allocator.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+radiative_env_release(ref_T* ref)
+{
+ struct sdis_radiative_env* radenv = NULL;
+ struct sdis_device* dev = NULL;
+ ASSERT(ref);
+ radenv = CONTAINER_OF(ref, struct sdis_radiative_env, ref);
+ dev = radenv->dev;
+ if(radenv->data) SDIS(data_ref_put(radenv->data));
+ MEM_RM(dev->allocator, radenv);
+ SDIS(device_ref_put(dev));
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+sdis_radiative_env_create
+ (struct sdis_device* dev,
+ const struct sdis_radiative_env_shader* shader,
+ struct sdis_data* data, /* Data sent to the shader. May be NULL */
+ struct sdis_radiative_env** out_radenv)
+{
+ struct sdis_radiative_env* radenv = NULL;
+ res_T res = RES_OK;
+
+ if(!dev || !shader || !out_radenv) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ radenv = MEM_CALLOC(dev->allocator, 1, sizeof(*radenv));
+ if(!radenv) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&radenv->ref);
+ SDIS(device_ref_get(dev));
+ radenv->dev = dev;
+ radenv->shader = *shader;
+ if(data) {
+ SDIS(data_ref_get(data));
+ radenv->data = data;
+ }
+
+exit:
+ if(out_radenv) *out_radenv = radenv;
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+sdis_radiative_env_ref_get(struct sdis_radiative_env* radenv)
+{
+ if(!radenv) return RES_BAD_ARG;
+ ref_get(&radenv->ref);
+ return RES_OK;
+}
+
+res_T
+sdis_radiative_env_ref_put(struct sdis_radiative_env* radenv)
+{
+ if(!radenv) return RES_BAD_ARG;
+ ref_put(&radenv->ref, radiative_env_release);
+ return RES_OK;
+}
+
+res_T
+sdis_radiative_env_get_shader
+ (struct sdis_radiative_env* radenv,
+ struct sdis_radiative_env_shader* shader)
+{
+ if(!radenv || !shader) return RES_BAD_ARG;
+ *shader = radenv->shader;
+ return RES_OK;
+}
+
+struct sdis_data*
+sdis_radiative_env_get_data(struct sdis_radiative_env* radenv)
+{
+ ASSERT(radenv);
+ return radenv->data;
+}
diff --git a/src/sdis_radiative_env_c.h b/src/sdis_radiative_env_c.h
@@ -0,0 +1,52 @@
+/* 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/>. */
+
+#ifndef SDIS_RADIATIVE_ENV_C_H
+#define SDIS_RADIATIVE_ENV_C_H
+
+#include "sdis.h"
+#include <rsys/ref_count.h>
+
+struct sdis_radiative_env {
+ struct sdis_radiative_env_shader shader;
+ struct sdis_data* data;
+
+ ref_T ref;
+ struct sdis_device* dev;
+};
+
+static INLINE double
+radiative_env_get_temperature
+ (const struct sdis_radiative_env* radenv,
+ const struct sdis_radiative_ray* ray)
+{
+ ASSERT(ray && d3_is_normalized(ray->dir));
+ return radenv && radenv->shader.temperature
+ ? radenv->shader.temperature(ray, radenv->data)
+ : SDIS_TEMPERATURE_NONE;
+}
+
+static INLINE double
+radiative_env_get_reference_temperature
+ (const struct sdis_radiative_env* radenv,
+ const struct sdis_radiative_ray* ray)
+{
+ ASSERT(ray && d3_is_normalized(ray->dir));
+ return radenv && radenv->shader.reference_temperature
+ ? radenv->shader.reference_temperature(ray, radenv->data)
+ : SDIS_TEMPERATURE_NONE;
+}
+
+#endif /* SDIS_RADIATIVE_ENV_C_H */
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -108,6 +108,7 @@ scene_release(ref_T * ref)
if(scn->senc2d_scn) SENC2D(scene_ref_put(scn->senc2d_scn));
if(scn->senc3d_scn) SENC3D(scene_ref_put(scn->senc3d_scn));
if(scn->source) SDIS(source_ref_put(scn->source));
+ if(scn->radenv) SDIS(radiative_env_ref_put(scn->radenv));
MEM_RM(dev->allocator, scn);
SDIS(device_ref_put(dev));
}
@@ -194,26 +195,6 @@ sdis_scene_set_fp_to_meter
}
res_T
-sdis_scene_get_ambient_radiative_temperature
- (const struct sdis_scene* scn,
- struct sdis_ambient_radiative_temperature* trad)
-{
- if(!scn || !trad) return RES_BAD_ARG;
- *trad = scn->trad;
- return RES_OK;
-}
-
-res_T
-sdis_scene_set_ambient_radiative_temperature
- (struct sdis_scene* scn,
- const struct sdis_ambient_radiative_temperature* trad)
-{
- if(!scn) return RES_BAD_ARG;
- scn->trad = *trad;
- return RES_OK;
-}
-
-res_T
sdis_scene_get_temperature_range
(const struct sdis_scene* scn,
double t_range[2])
@@ -425,6 +406,16 @@ sdis_scene_get_source(struct sdis_scene* scn, struct sdis_source** source)
return RES_OK;
}
+res_T
+sdis_scene_get_radiative_env
+ (struct sdis_scene* scn,
+ struct sdis_radiative_env** radenv)
+{
+ if(!scn || !radenv) return RES_BAD_ARG;
+ *radenv = scn->radenv;
+ return RES_OK;
+}
+
/*******************************************************************************
* Local miscellaneous function
******************************************************************************/
@@ -463,6 +454,7 @@ scene_compute_hash(const struct sdis_scene* scn, hash256_T hash)
{
struct sha256_ctx sha256_ctx;
size_t iprim, nprims;
+ int has_radenv = 0;
res_T res = RES_OK;
ASSERT(scn && hash);
@@ -476,7 +468,9 @@ scene_compute_hash(const struct sdis_scene* scn, hash256_T hash)
#define SHA256_UPD(Var, Nb) \
sha256_ctx_update(&sha256_ctx, (const char*)(Var), sizeof(*Var)*(Nb))
- SHA256_UPD(&scn->trad.reference, 1);
+ has_radenv = scn->radenv != NULL;
+
+ SHA256_UPD(&has_radenv, 1);
SHA256_UPD(&scn->tmax, 1);
SHA256_UPD(&scn->fp_to_meter, 1);
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -965,7 +965,6 @@ XD(scene_create)
SDIS(device_ref_get(dev));
scn->dev = dev;
scn->fp_to_meter = args->fp_to_meter;
- scn->trad = args->trad;
scn->tmin = args->t_range[0];
scn->tmax = args->t_range[1];
scn->outer_enclosure_id = UINT_MAX;
@@ -980,6 +979,11 @@ XD(scene_create)
scn->source = args->source;
}
+ if(args->radenv) {
+ SDIS(radiative_env_ref_get(args->radenv));
+ scn->radenv = args->radenv;
+ }
+
res = XD(run_analyze)
(scn,
args->nprimitives,
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -221,11 +221,11 @@ struct sdis_scene {
unsigned outer_enclosure_id;
double fp_to_meter;
- struct sdis_ambient_radiative_temperature trad;
double tmin; /* Minimum temperature of the system (In Kelvin) */
double tmax; /* Maximum temperature of the system (In Kelvin) */
struct sdis_source* source; /* External source. May be NULL */
+ struct sdis_radiative_env* radenv; /* Radiative environment. May be NULL */
ref_T ref;
struct sdis_device* dev;
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -784,7 +784,7 @@ XD(solve_boundary_flux)
if(res_local!= RES_OK) { ATOMIC_SET(&res, res_local); continue; }
/* Fetch interface parameters */
- epsilon = interface_side_get_emissivity(interf, &frag);
+ epsilon = interface_side_get_emissivity(interf, SDIS_INTERN_SOURCE_ID, &frag);
hc = interface_get_convection_coef(interf, &frag);
Tref = interface_side_get_reference_temperature(interf, &frag);
if(epsilon <= 0) {
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -880,7 +880,8 @@ XD(solve_probe_boundary_flux)
/* Compute hr and hc */
frag_local.time = time;
frag_local.side = fluid_side;
- epsilon = interface_side_get_emissivity(interf, &frag_local);
+ epsilon = interface_side_get_emissivity
+ (interf, SDIS_INTERN_SOURCE_ID, &frag_local);
Tref = interface_side_get_reference_temperature(interf, &frag_local);
hc = interface_get_convection_coef(interf, &frag_local);
if(epsilon <= 0) {
diff --git a/src/sdis_source.c b/src/sdis_source.c
@@ -18,6 +18,7 @@
#include "sdis_log.h"
#include "sdis_source_c.h"
+#include <rsys/free_list.h>
#include <rsys/mem_allocator.h>
#include <rsys/ref_count.h>
@@ -26,6 +27,7 @@
struct sdis_source {
struct sdis_spherical_source_create_args spherical;
+ struct fid id; /* Unique identifier of the source */
struct sdis_device* dev;
ref_T ref;
};
@@ -69,6 +71,7 @@ release_source(ref_T* ref)
ASSERT(ref);
dev = src->dev;
if(src->spherical.data) SDIS(data_ref_put(src->spherical.data));
+ flist_name_del(&dev->source_names, src->id);
MEM_RM(dev->allocator, src);
SDIS(device_ref_put(dev));
}
@@ -100,6 +103,8 @@ sdis_spherical_source_create
if(args->data) SDIS(data_ref_get(args->data));
src->spherical = *args;
src->dev = dev;
+ src->id = flist_name_add(&dev->source_names);
+ flist_name_get(&dev->source_names, src->id)->mem = src;
exit:
if(out_src) *out_src = src;
@@ -132,6 +137,13 @@ sdis_source_get_power(struct sdis_source* src, const double time /* [s] */)
return source_get_power(src, time);
}
+unsigned
+sdis_source_get_id(const struct sdis_source* source)
+{
+ ASSERT(source);
+ return source->id.index;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -212,16 +212,22 @@ interface_get_convection_coef
static double
interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(data != NULL && frag != NULL);
return ((const struct interfac*)sdis_data_cget(data))->emissivity;
}
static double
interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(data != NULL && frag != NULL);
return ((const struct interfac*)sdis_data_cget(data))->specular_fraction;
}
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -192,10 +192,13 @@ interface_get_convection_coef
static double
interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
const struct interfac* interf;
double e = -1;
+ (void)source_id;
CHK(data != NULL && frag != NULL);
interf = sdis_data_cget(data);
switch(frag->side) {
@@ -208,10 +211,13 @@ interface_get_emissivity
static double
interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
const struct interfac* interf;
double f = -1;
+ (void)source_id;
CHK(data != NULL && frag != NULL);
interf = sdis_data_cget(data);
switch(frag->side) {
diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c
@@ -122,16 +122,22 @@ interface_get_convection_coef
static double
interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(frag && data);
return 0;
}
static double
interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(frag && data);
return 0;
}
diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c
@@ -131,16 +131,22 @@ interface_get_convection_coef
static double
interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(frag && data);
return 0;
}
static double
interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(frag && data);
return 0;
}
diff --git a/src/test_sdis_draw_external_flux.c b/src/test_sdis_draw_external_flux.c
@@ -144,9 +144,11 @@ interface_get_temperature
static double
interface_get_emissivity
(const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
struct sdis_data* data)
{
- (void)frag, (void)data;/* Avoid the "unused variable" warning */
+ /* Avoid the "unused variable" warning */
+ (void)frag, (void)source_id, (void)data;
return EMISSIVITY;
}
@@ -204,6 +206,39 @@ create_source(struct sdis_device* sdis)
}
/*******************************************************************************
+ * Radiative environment
+ ******************************************************************************/
+static double
+radenv_get_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray, (void)data;
+ return T_RAD;
+}
+
+static double
+radenv_get_reference_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray, (void)data;
+ return T_REF;
+}
+
+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;
+}
+
+/*******************************************************************************
* Scene, i.e. the system to simulate
******************************************************************************/
struct scene_context {
@@ -245,7 +280,8 @@ create_scene
(struct sdis_device* sdis,
const struct mesh* mesh,
struct sdis_interface* interf,
- struct sdis_source* source)
+ 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;
@@ -259,11 +295,10 @@ create_scene
scn_args.get_position = scene_get_position;
scn_args.nprimitives = mesh_ntriangles(mesh);
scn_args.nvertices = mesh_nvertices(mesh);
- scn_args.trad.temperature = T_RAD;
- scn_args.trad.reference = T_REF;
scn_args.t_range[0] = MMIN(T_RAD, T_REF);
scn_args.t_range[1] = MMAX(T_RAD, T_REF);
scn_args.source = source;
+ scn_args.radenv = radenv;
scn_args.context = &context;
OK(sdis_scene_create(sdis, &scn_args, &scn));
return scn;
@@ -342,6 +377,7 @@ main(int argc, char** argv)
struct sdis_medium* solid = NULL;
struct sdis_scene* scn = NULL;
struct sdis_source* source = NULL;
+ struct sdis_radiative_env* radenv = NULL;
/* Miscellaneous */
struct mesh mesh = MESH_NULL;
@@ -357,7 +393,8 @@ main(int argc, char** argv)
fluid = create_fluid(dev);
interf = create_interface(dev, fluid, solid);
source = create_source(dev);
- scn = create_scene(dev, &mesh, interf, source);
+ radenv = create_radenv(dev);
+ scn = create_scene(dev, &mesh, interf, source, radenv);
cam = create_camera(dev);
draw(scn, cam);
@@ -375,6 +412,7 @@ main(int argc, char** argv)
OK(sdis_medium_ref_put(solid));
OK(sdis_scene_ref_put(scn));
OK(sdis_source_ref_put(source));
+ OK(sdis_radiative_env_ref_put(radenv));
CHK(mem_allocated_size() == 0);
return 0;
}
diff --git a/src/test_sdis_external_flux.c b/src/test_sdis_external_flux.c
@@ -45,7 +45,7 @@
*/
#define T_FLUID 300.0 /* [K] */
-#define T_REF 0 /* [K] */
+#define T_REF 300.0 /* [K] */
/*******************************************************************************
* Geometries
@@ -196,18 +196,29 @@ INTERF_PROP(reference_temperature, T_REF) /* [K] */
static double \
interface_get_##Prop \
(const struct sdis_interface_fragment* frag, \
+ const unsigned source_id, \
struct sdis_data* data) \
{ \
struct interface* interf_data = NULL; \
- (void)frag; /* Avoid the "unused variable" warning */ \
+ (void)frag, (void)source_id; /* Avoid the "unused variable" warning */ \
interf_data = sdis_data_get(data); \
- return interf_data->Prop; \
+ return source_id == SDIS_INTERN_SOURCE_ID ? 0 : interf_data->Prop; \
}
INTERF_PROP(emissivity)
INTERF_PROP(specular_fraction)
-INTERF_PROP(convection_coef) /* [W/m^2/K] */
#undef INTERF_PROP
+static double /* [W/m^2/K] */
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ struct interface* interf_data = NULL;
+ (void)frag; /* Avoid the "unused variable" warning */
+ interf_data = sdis_data_get(data);
+ return interf_data->convection_coef;
+}
+
static struct sdis_interface*
create_interface
(struct sdis_device* sdis,
@@ -283,6 +294,39 @@ create_source(struct sdis_device* sdis)
}
/*******************************************************************************
+ * Radiative environment
+ ******************************************************************************/
+static double
+radenv_get_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray, (void)data;
+ return 0; /* [K] */
+}
+
+static double
+radenv_get_reference_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray, (void)data;
+ return T_REF; /* [K] */
+}
+
+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;
+}
+
+/*******************************************************************************
* Scene
******************************************************************************/
struct scene_context {
@@ -350,7 +394,8 @@ create_scene_2d
(struct sdis_device* sdis,
struct sdis_interface* interf_ground,
struct sdis_interface* interf_wall,
- struct sdis_source* source)
+ struct sdis_source* source,
+ struct sdis_radiative_env* radenv)
{
struct sdis_scene* scn = NULL;
struct sdis_source* src = NULL;
@@ -365,11 +410,10 @@ create_scene_2d
scn_args.get_position = scene_get_position_2d;
scn_args.nprimitives = nsegments;
scn_args.nvertices = nvertices_2d;
- scn_args.trad.temperature = 0; /* [K] */
- scn_args.trad.reference = T_REF; /* [K] */
- scn_args.t_range[0] = 0; /* [K] */
- scn_args.t_range[1] = 0; /* [K] */
+ scn_args.t_range[0] = T_REF; /* [K] */
+ scn_args.t_range[1] = T_REF; /* [K] */
scn_args.source = source;
+ scn_args.radenv = radenv;
scn_args.context = &context;
OK(sdis_scene_2d_create(sdis, &scn_args, &scn));
@@ -386,7 +430,8 @@ create_scene_3d
(struct sdis_device* sdis,
struct sdis_interface* interf_ground,
struct sdis_interface* interf_wall,
- struct sdis_source* source)
+ struct sdis_source* source,
+ struct sdis_radiative_env* radenv)
{
struct sdis_scene* scn = NULL;
struct sdis_source* src = NULL;
@@ -401,11 +446,10 @@ create_scene_3d
scn_args.get_position = scene_get_position_3d;
scn_args.nprimitives = ntriangles;
scn_args.nvertices = nvertices_3d;
- scn_args.trad.temperature = 0; /* [K] */
- scn_args.trad.reference = T_REF; /* [K] */
scn_args.t_range[0] = 0; /* [K] */
scn_args.t_range[1] = 0; /* [K] */
scn_args.source = source;
+ scn_args.radenv = radenv;
scn_args.context = &context;
OK(sdis_scene_create(sdis, &scn_args, &scn));
@@ -491,9 +535,10 @@ main(int argc, char** argv)
struct sdis_medium* solid = NULL;
struct sdis_interface* interf_ground = NULL;
struct sdis_interface* interf_wall = NULL;
- struct sdis_source* src = NULL;
+ struct sdis_radiative_env* radenv = NULL;
struct sdis_scene* scn_2d = NULL;
struct sdis_scene* scn_3d = NULL;
+ struct sdis_source* src = NULL;
struct interface* ground_interf_data = NULL;
int is_master_process = 0;
@@ -508,8 +553,9 @@ main(int argc, char** argv)
interf_wall = create_interface
(dev, solid, fluid, 1/*emissivity*/, 10/*h*/, NULL);
src = create_source(dev);
- scn_2d = create_scene_2d(dev, interf_ground, interf_wall, src);
- scn_3d = create_scene_3d(dev, interf_ground, interf_wall, src);
+ radenv = create_radenv(dev);
+ scn_2d = create_scene_2d(dev, interf_ground, interf_wall, src, radenv);
+ scn_3d = create_scene_3d(dev, interf_ground, interf_wall, src, radenv);
ground_interf_data->specular_fraction = 0; /* Lambertian */
check(scn_2d, 10000, 375.88, is_master_process);
@@ -523,6 +569,7 @@ main(int argc, char** argv)
OK(sdis_medium_ref_put(solid));
OK(sdis_interface_ref_put(interf_ground));
OK(sdis_interface_ref_put(interf_wall));
+ OK(sdis_radiative_env_ref_put(radenv));
OK(sdis_source_ref_put(src));
OK(sdis_scene_ref_put(scn_2d));
OK(sdis_scene_ref_put(scn_3d));
diff --git a/src/test_sdis_flux2.c b/src/test_sdis_flux2.c
@@ -281,17 +281,23 @@ interface_get_convection_coef
static double
interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
const struct interf* interf = sdis_data_cget(data);
+ (void)source_id;
CHK(frag && interf);
return interf->emissivity;
}
static double
interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(frag && data);
return 0; /* Unused */
}
@@ -347,19 +353,53 @@ create_interface
}
/*******************************************************************************
+ * Create the radiative environment
+ ******************************************************************************/
+static double
+radenv_get_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray, (void)data;
+ return 320; /* [K] */
+}
+
+static double
+radenv_get_reference_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray, (void)data;
+ return 300; /* [K] */
+}
+
+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;
+}
+
+/*******************************************************************************
* Create scene
******************************************************************************/
static void
create_scene_3d
(struct sdis_device* dev,
struct sdis_interface* interfaces[INTERFACES_COUNT__],
+ struct sdis_radiative_env* radenv,
struct sdis_scene** scn)
{
struct geometry geom;
struct sdis_interface* prim_interfaces[32];
struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
- CHK(dev && interfaces && scn);
+ CHK(dev && interfaces && radenv && scn);
/* Setup the per primitive interface of the solid medium */
prim_interfaces[0] = prim_interfaces[1] = interfaces[ADIABATIC];
@@ -388,8 +428,7 @@ create_scene_3d
scn_args.t_range[0] = 300;
scn_args.t_range[1] = 300;
scn_args.context = &geom;
- scn_args.trad.temperature = 320;
- scn_args.trad.reference = 300;
+ scn_args.radenv = radenv;
OK(sdis_scene_create(dev, &scn_args, scn));
}
@@ -430,6 +469,7 @@ int
main(int argc, char** argv)
{
struct sdis_device* dev = NULL;
+ struct sdis_radiative_env* radenv = NULL;
struct sdis_scene* scn_3d = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* dummy = NULL;
@@ -465,6 +505,8 @@ main(int argc, char** argv)
OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev));
+ radenv = create_radenv(dev);
+
/* Solid medium */
solid_props.lambda = 1.15;
solid_props.rho = 1700;
@@ -517,13 +559,14 @@ main(int argc, char** argv)
create_interface
(dev, solid, fluid2, &interf_props, &interfaces[SOLID_FLUID]);
- create_scene_3d(dev, interfaces, &scn_3d);
+ create_scene_3d(dev, interfaces, radenv, &scn_3d);
FOR_EACH(iprobe, 0, nprobes) {
check(scn_3d, &probes[iprobe]);
}
/* Release memory */
+ OK(sdis_radiative_env_ref_put(radenv));
OK(sdis_scene_ref_put(scn_3d));
OK(sdis_medium_ref_put(solid));
OK(sdis_medium_ref_put(dummy));
diff --git a/src/test_sdis_interface.c b/src/test_sdis_interface.c
@@ -88,18 +88,18 @@ main(int argc, char** argv)
OK(CREATE(dev, solid, solid, &shader, NULL, &interf));
OK(sdis_interface_ref_put(interf));
- shader.back.emissivity = dummy_interface_getter;
+ shader.back.emissivity = dummy_radiative_interface_getter;
OK(CREATE(dev, solid, fluid, &shader, NULL, &interf));
OK(sdis_interface_ref_put(interf));
- shader.back.specular_fraction = dummy_interface_getter;
+ shader.back.specular_fraction = dummy_radiative_interface_getter;
OK(CREATE(dev, solid, fluid, &shader, NULL, &interf));
OK(sdis_interface_ref_put(interf));
shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
- shader.front.emissivity = dummy_interface_getter;
+ shader.front.emissivity = dummy_radiative_interface_getter;
OK(CREATE(dev, solid, fluid, &shader, NULL, &interf)); /* Warning */
OK(sdis_interface_ref_put(interf));
shader.front.emissivity = NULL;
- shader.front.specular_fraction = dummy_interface_getter;
+ shader.front.specular_fraction = dummy_radiative_interface_getter;
OK(CREATE(dev, solid, fluid, &shader, NULL, &interf)); /* Warning */
OK(sdis_interface_ref_put(interf));
shader.front.specular_fraction = NULL;
diff --git a/src/test_sdis_picard.c b/src/test_sdis_picard.c
@@ -306,18 +306,24 @@ interface_get_convection_coef
static double
interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
const struct interf* interf = sdis_data_cget(data);
+ (void)source_id;
CHK(frag && interf);
return interf->emissivity;
}
static double
interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
const struct interf* interf = sdis_data_cget(data);
+ (void)source_id;
CHK(frag && interf);
return interf->specular_fraction;
}
@@ -370,6 +376,52 @@ create_interface
}
/*******************************************************************************
+ * Radiative environment
+ ******************************************************************************/
+struct radenv {
+ double temperature; /* [K] */
+ double reference; /* [K] */
+};
+
+static double
+radenv_get_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray;
+ return ((const struct radenv*)sdis_data_cget(data))->temperature;
+}
+
+static double
+radenv_get_reference_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray;
+ return ((const struct radenv*)sdis_data_cget(data))->reference;
+}
+
+static struct sdis_radiative_env*
+create_radenv(struct sdis_device* dev)
+{
+ struct sdis_radiative_env_shader shader = SDIS_RADIATIVE_ENV_SHADER_NULL;
+ struct sdis_radiative_env* radenv = NULL;
+ struct sdis_data* data = NULL;
+ struct radenv* env = NULL;
+
+ OK(sdis_data_create(dev, sizeof(struct radenv), ALIGNOF(radenv), NULL, &data));
+ env = sdis_data_get(data);
+ env->temperature = 300;
+ env->reference = 300;
+
+ shader.temperature = radenv_get_temperature;
+ shader.reference_temperature = radenv_get_reference_temperature;
+ OK(sdis_radiative_env_create(dev, &shader, data, &radenv));
+ OK(sdis_data_ref_put(data));
+ return radenv;
+}
+
+/*******************************************************************************
* Helper functions
******************************************************************************/
struct reference_result {
@@ -482,13 +534,14 @@ static void
create_scene_3d
(struct sdis_device* dev,
struct sdis_interface* interfaces[INTERFACES_COUNT__],
+ struct sdis_radiative_env* radenv,
struct sdis_scene** scn)
{
struct geometry geom;
struct sdis_interface* prim_interfaces[32];
struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
- CHK(dev && interfaces && scn);
+ CHK(dev && interfaces && radenv && scn);
/* Setup the per primitive interface of the solid medium */
prim_interfaces[0] = prim_interfaces[1] = interfaces[ADIABATIC];
@@ -516,6 +569,7 @@ create_scene_3d
scn_args.nvertices = nvertices_3d;
scn_args.t_range[0] = 280;
scn_args.t_range[1] = 350;
+ scn_args.radenv = radenv;
scn_args.context = &geom;
OK(sdis_scene_create(dev, &scn_args, scn));
}
@@ -524,13 +578,14 @@ static void
create_scene_2d
(struct sdis_device* dev,
struct sdis_interface* interfaces[INTERFACES_COUNT__],
+ struct sdis_radiative_env* radenv,
struct sdis_scene** scn)
{
struct geometry geom;
struct sdis_interface* prim_interfaces[10/*#segment*/];
struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
- CHK(dev && interfaces && scn);
+ CHK(dev && interfaces && radenv && scn);
/* Setup the per primitive interface of the solid medium */
prim_interfaces[0] = interfaces[ADIABATIC];
@@ -554,6 +609,7 @@ create_scene_2d
scn_args.nvertices = nvertices_2d;
scn_args.t_range[0] = 280;
scn_args.t_range[1] = 350;
+ scn_args.radenv = radenv;
scn_args.context = &geom;
OK(sdis_scene_2d_create(dev, &scn_args, scn));
}
@@ -567,14 +623,15 @@ main(int argc, char** argv)
FILE* stream = NULL;
struct sdis_device* dev = NULL;
+ struct sdis_radiative_env* radenv = NULL;
struct sdis_scene* scn_2d = NULL;
struct sdis_scene* scn_3d = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid = NULL;
struct sdis_medium* dummy = NULL;
struct sdis_interface* interfaces[INTERFACES_COUNT__];
- struct sdis_ambient_radiative_temperature amb_rad_temp;
+ struct radenv* radenv_props = NULL;
struct solid solid_props;
struct solid* psolid_props;
struct reference_result ref = REFERENCE_RESULT_NULL;
@@ -588,6 +645,9 @@ main(int argc, char** argv)
OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev));
+ radenv = create_radenv(dev);
+ radenv_props = sdis_data_get(sdis_radiative_env_get_data(radenv));
+
/* Solid medium */
solid_props.lambda = 1.15;
solid_props.rho = 1000;
@@ -637,8 +697,8 @@ main(int argc, char** argv)
pinterf_props[i] = sdis_data_get(sdis_interface_get_data(interfaces[i]));
}
- create_scene_2d(dev, interfaces, &scn_2d);
- create_scene_3d(dev, interfaces, &scn_3d);
+ create_scene_2d(dev, interfaces, radenv, &scn_2d);
+ create_scene_3d(dev, interfaces, radenv, &scn_3d);
CHK((stream = tmpfile()) != NULL);
@@ -650,10 +710,8 @@ main(int argc, char** argv)
pinterf_props[SOLID_FLUID_mX]->Tref = 300;
pinterf_props[SOLID_FLUID_pX]->Tref = 300;
pinterf_props[BOUNDARY_pX]->Tref = 300;
- amb_rad_temp.temperature = 280;
- amb_rad_temp.reference = 300;
- OK(sdis_scene_set_ambient_radiative_temperature(scn_2d, &amb_rad_temp));
- OK(sdis_scene_set_ambient_radiative_temperature(scn_3d, &amb_rad_temp));
+ radenv_props->temperature = 280;
+ radenv_props->reference = 300;
test_picard(scn_2d, 1/*Picard order*/, &ref);
test_picard(scn_3d, 1/*Picard order*/, &ref);
printf("\n");
@@ -666,10 +724,8 @@ main(int argc, char** argv)
pinterf_props[SOLID_FLUID_mX]->Tref = ref.T1;
pinterf_props[SOLID_FLUID_pX]->Tref = ref.T2;
pinterf_props[BOUNDARY_pX]->Tref = 350;
- amb_rad_temp.temperature = 280;
- amb_rad_temp.reference = 280;
- OK(sdis_scene_set_ambient_radiative_temperature(scn_2d, &amb_rad_temp));
- OK(sdis_scene_set_ambient_radiative_temperature(scn_3d, &amb_rad_temp));
+ radenv_props->temperature = 280;
+ radenv_props->reference = 280;
test_picard(scn_2d, 1/*Picard order*/, &ref);
test_picard(scn_3d, 1/*Picard order*/, &ref);
printf("\n");
@@ -682,10 +738,8 @@ main(int argc, char** argv)
pinterf_props[SOLID_FLUID_mX]->Tref = 300;
pinterf_props[SOLID_FLUID_pX]->Tref = 300;
pinterf_props[BOUNDARY_pX]->Tref = 300;
- amb_rad_temp.temperature = 280;
- amb_rad_temp.reference = 300;
- OK(sdis_scene_set_ambient_radiative_temperature(scn_2d, &amb_rad_temp));
- OK(sdis_scene_set_ambient_radiative_temperature(scn_3d, &amb_rad_temp));
+ radenv_props->temperature = 280;
+ radenv_props->reference = 300;
test_picard(scn_2d, 2/*Picard order*/, &ref);
test_picard(scn_3d, 2/*Picard order*/, &ref);
printf("\n");
@@ -705,10 +759,8 @@ main(int argc, char** argv)
pinterf_props[SOLID_FLUID_mX]->Tref = 350;
pinterf_props[SOLID_FLUID_pX]->Tref = 450;
pinterf_props[BOUNDARY_pX]->Tref = pinterf_props[BOUNDARY_pX]->temperature;
- amb_rad_temp.temperature = t_range[0];
- amb_rad_temp.reference = t_range[0];
- OK(sdis_scene_set_ambient_radiative_temperature(scn_2d, &amb_rad_temp));
- OK(sdis_scene_set_ambient_radiative_temperature(scn_3d, &amb_rad_temp));
+ radenv_props->temperature = t_range[0];
+ radenv_props->reference = t_range[0];
test_picard(scn_2d, 3/*Picard order*/, &ref);
test_picard(scn_3d, 3/*Picard order*/, &ref);
register_heat_paths(scn_2d, 3/*Picard order*/, stream);
@@ -733,10 +785,8 @@ main(int argc, char** argv)
pinterf_props[SOLID_FLUID_mX]->Tref = 300;
pinterf_props[SOLID_FLUID_pX]->Tref = 300;
pinterf_props[BOUNDARY_pX]->Tref = 300;
- amb_rad_temp.temperature = t_range[0];
- amb_rad_temp.reference = 300;
- OK(sdis_scene_set_ambient_radiative_temperature(scn_2d, &amb_rad_temp));
- OK(sdis_scene_set_ambient_radiative_temperature(scn_3d, &amb_rad_temp));
+ radenv_props->temperature = t_range[0];
+ radenv_props->reference = 300;
test_picard(scn_2d, 1/*Picard order*/, &ref);
test_picard(scn_3d, 1/*Picard order*/, &ref);
printf("\n");
@@ -749,15 +799,14 @@ main(int argc, char** argv)
pinterf_props[SOLID_FLUID_mX]->Tref = ref.T1;
pinterf_props[SOLID_FLUID_pX]->Tref = ref.T2;
pinterf_props[BOUNDARY_pX]->Tref = 350;
- amb_rad_temp.temperature = 280;
- amb_rad_temp.reference = 280;
- OK(sdis_scene_set_ambient_radiative_temperature(scn_2d, &amb_rad_temp));
- OK(sdis_scene_set_ambient_radiative_temperature(scn_3d, &amb_rad_temp));
+ radenv_props->temperature = 280;
+ radenv_props->reference = 280;
test_picard(scn_2d, 1/*Picard order*/, &ref);
test_picard(scn_3d, 1/*Picard order*/, &ref);
printf("\n");
/* Release memory */
+ OK(sdis_radiative_env_ref_put(radenv));
OK(sdis_scene_ref_put(scn_2d));
OK(sdis_scene_ref_put(scn_3d));
OK(sdis_interface_ref_put(interfaces[ADIABATIC]));
diff --git a/src/test_sdis_radiative_env.c b/src/test_sdis_radiative_env.c
@@ -0,0 +1,107 @@
+/* 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 "sdis.h"
+#include "test_sdis_utils.h"
+
+#include <rsys/rsys.h>
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static double
+radenv_get_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray, (void)data;
+ return 300;
+}
+
+static double
+radenv_get_reference_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray, (void)data;
+ return 300;
+}
+
+static void
+check_api(struct sdis_device* sdis)
+{
+ struct sdis_radiative_env_shader shader = SDIS_RADIATIVE_ENV_SHADER_NULL;
+ struct sdis_radiative_env_shader shader2 = SDIS_RADIATIVE_ENV_SHADER_NULL;
+ struct sdis_radiative_env* radenv = NULL;
+ struct sdis_data* data = NULL;
+
+ CHK(sdis != NULL);
+ BA(sdis_radiative_env_create(NULL, &shader, NULL, &radenv));
+ BA(sdis_radiative_env_create(sdis, NULL, NULL, &radenv));
+ OK(sdis_radiative_env_create(sdis, &shader, NULL, &radenv));
+
+ BA(sdis_radiative_env_get_shader(NULL, &shader2));
+ BA(sdis_radiative_env_get_shader(radenv, NULL));
+ OK(sdis_radiative_env_get_shader(radenv, &shader2));
+
+ CHK(shader2.temperature == shader.temperature);
+ CHK(shader2.reference_temperature == shader.reference_temperature);
+
+ CHK(sdis_radiative_env_get_data(radenv) == NULL);
+
+ BA(sdis_radiative_env_ref_get(NULL));
+ OK(sdis_radiative_env_ref_get(radenv));
+ BA(sdis_radiative_env_ref_put(NULL));
+ OK(sdis_radiative_env_ref_put(radenv));
+ OK(sdis_radiative_env_ref_put(radenv));
+
+ shader.temperature = radenv_get_temperature;
+ shader.reference_temperature = radenv_get_reference_temperature;
+
+ OK(sdis_data_create(sdis, sizeof(uint32_t), ALIGNOF(uint32_t), NULL, &data));
+ *((uint32_t*)sdis_data_get(data)) = 0xD3CAFBAD;
+
+ OK(sdis_radiative_env_create(sdis, &shader, data, &radenv));
+ OK(sdis_data_ref_put(data));
+
+ OK(sdis_radiative_env_get_shader(radenv, &shader2));
+
+ CHK(shader2.temperature == shader.temperature);
+ CHK(shader2.reference_temperature == shader.reference_temperature);
+
+ CHK(sdis_radiative_env_get_data(radenv) == data);
+ data = sdis_radiative_env_get_data(radenv);
+ CHK(*((uint32_t*)sdis_data_get(data)) == 0xD3CAFBAD);
+
+ OK(sdis_radiative_env_ref_put(radenv));
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct sdis_device* sdis = NULL;
+ (void)argc, (void)argv;
+
+ OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &sdis));
+
+ check_api(sdis);
+
+ OK(sdis_device_ref_put(sdis));
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c
@@ -87,7 +87,10 @@ get_interface(const size_t itri, struct sdis_interface** bound, void* context)
}
static void
-test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
+test_scene_3d
+ (struct sdis_device* dev,
+ struct sdis_interface* interf,
+ struct sdis_radiative_env* in_radenv)
{
size_t duplicated_indices[] = { 0, 1, 2, 0, 2, 1 };
size_t degenerated_indices[] = { 0, 1, 1 };
@@ -102,6 +105,7 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL;
struct sdis_scene* scn = NULL;
struct sdis_device* dev2 = NULL;
+ struct sdis_radiative_env* radenv = NULL;
size_t ntris, npos;
size_t iprim;
size_t i;
@@ -266,15 +270,31 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
/* No 2D available */
BA(sdis_scene_get_senc2d_scene(scn, &scn2d));
+ BA(sdis_scene_get_radiative_env(NULL, &radenv));
+ BA(sdis_scene_get_radiative_env(scn, NULL));
+ OK(sdis_scene_get_radiative_env(scn, &radenv));
+ CHK(radenv == NULL);
+
BA(sdis_scene_ref_get(NULL));
OK(sdis_scene_ref_get(scn));
BA(sdis_scene_ref_put(NULL));
OK(sdis_scene_ref_put(scn));
OK(sdis_scene_ref_put(scn));
+
+ scn_args.radenv = in_radenv;
+ OK(sdis_scene_create(dev, &scn_args, &scn));
+ OK(sdis_scene_get_radiative_env(scn, &radenv));
+ CHK(radenv == in_radenv);
+ CHK(radenv != NULL);
+
+ OK(sdis_scene_ref_put(scn));
}
static void
-test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
+test_scene_2d
+ (struct sdis_device* dev,
+ struct sdis_interface* interf,
+ struct sdis_radiative_env* in_radenv)
{
size_t duplicated_indices[] = { 0, 1, 1, 0 };
size_t degenerated_indices[] = { 0, 0 };
@@ -283,8 +303,7 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
struct sdis_scene_find_closest_point_args closest_pt_args =
SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL;
- struct sdis_ambient_radiative_temperature trad =
- SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
+ struct sdis_radiative_env* radenv = NULL;
double lower[2], upper[2];
double t_range[2];
double u0, u1, u2, pos[2], pos1[2];
@@ -379,30 +398,6 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
OK(sdis_scene_get_fp_to_meter(scn, &fp));
CHK(fp == 2);
- BA(sdis_scene_get_ambient_radiative_temperature(NULL, NULL));
- 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));
- 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;
- BA(sdis_scene_set_ambient_radiative_temperature(NULL, &trad));
- OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
- trad = SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
- OK(sdis_scene_get_ambient_radiative_temperature(scn, &trad));
- CHK(trad.temperature == 100);
- CHK(trad.reference == 110);
-
BA(sdis_scene_get_temperature_range(NULL, NULL));
BA(sdis_scene_get_temperature_range(scn, NULL));
BA(sdis_scene_get_temperature_range(NULL, t_range));
@@ -511,11 +506,25 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
/* No 3D available */
BA(sdis_scene_get_senc3d_scene(scn, &scn3d));
+ BA(sdis_scene_get_radiative_env(NULL, NULL));
+ BA(sdis_scene_get_radiative_env(scn, NULL));
+ BA(sdis_scene_get_radiative_env(NULL, &radenv));
+ OK(sdis_scene_get_radiative_env(scn, &radenv));
+ CHK(radenv == NULL);
+
BA(sdis_scene_ref_get(NULL));
OK(sdis_scene_ref_get(scn));
BA(sdis_scene_ref_put(NULL));
OK(sdis_scene_ref_put(scn));
OK(sdis_scene_ref_put(scn));
+
+ scn_args.radenv = in_radenv;
+ OK(sdis_scene_2d_create(dev, &scn_args, &scn));
+ OK(sdis_scene_get_radiative_env(scn, &radenv));
+ CHK(radenv == in_radenv);
+ CHK(radenv != NULL);
+
+ OK(sdis_scene_ref_put(scn));
}
int
@@ -525,9 +534,11 @@ main(int argc, char** argv)
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid = NULL;
struct sdis_interface* interf = NULL;
+ struct sdis_radiative_env* radenv = NULL;
struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_radiative_env_shader ray_shader = DUMMY_RAY_SHADER;
(void)argc, (void)argv;
interface_shader.convection_coef = DUMMY_INTERFACE_SHADER.convection_coef;
@@ -538,15 +549,17 @@ main(int argc, char** argv)
OK(sdis_solid_create(dev, &solid_shader, NULL, &solid));
OK(sdis_interface_create
(dev, solid, fluid, &interface_shader, NULL, &interf));
+ OK(sdis_radiative_env_create(dev, &ray_shader, NULL, &radenv));
OK(sdis_medium_ref_put(solid));
OK(sdis_medium_ref_put(fluid));
- test_scene_3d(dev, interf);
- test_scene_2d(dev, interf);
+ test_scene_3d(dev, interf, radenv);
+ test_scene_2d(dev, interf, radenv);
OK(sdis_device_ref_put(dev));
OK(sdis_interface_ref_put(interf));
+ OK(sdis_radiative_env_ref_put(radenv));
CHK(mem_allocated_size() == 0);
return 0;
diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c
@@ -156,9 +156,12 @@ interface_get_temperature
static double
interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
const struct interf* interf = sdis_data_cget(data);
+ (void)source_id;
CHK(frag && data);
return interf->emissivity;
}
@@ -182,6 +185,39 @@ interface_get_reference_temperature
}
/*******************************************************************************
+ * Radiative environment
+ ******************************************************************************/
+static double
+radenv_get_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray, (void)data;
+ return Trad;
+}
+
+static double
+radenv_get_reference_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray, (void)data;
+ return Trad;
+}
+
+static struct sdis_radiative_env*
+create_radenv(struct sdis_device* dev)
+{
+ 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(dev, &shader, NULL, &radenv));
+ return radenv;
+}
+
+/*******************************************************************************
* Helper function
******************************************************************************/
static void
@@ -235,6 +271,7 @@ main(int argc, char** argv)
struct sdis_interface* interf_adiabatic = NULL;
struct sdis_interface* interf_Tb = NULL;
struct sdis_interface* interf_H = NULL;
+ struct sdis_radiative_env* radenv = NULL;
struct sdis_scene* box_scn = NULL;
struct sdis_scene* square_scn = NULL;
struct sdis_estimator* estimator = NULL;
@@ -250,7 +287,7 @@ main(int argc, char** argv)
struct sdis_solve_boundary_flux_args bound_args =
SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT;
struct interf* interf_props = NULL;
- struct fluid* fluid_param;
+ struct fluid* fluid_args = NULL;
struct ssp_rng* rng = NULL;
enum sdis_estimator_type type;
double pos[3];
@@ -260,12 +297,13 @@ main(int argc, char** argv)
(void)argc, (void)argv;
create_default_device(&argc, &argv, &is_master_process, &dev);
+ radenv = create_radenv(dev);
/* Create the fluid medium */
OK(sdis_data_create
(dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
- fluid_param = sdis_data_get(data);
- fluid_param->temperature = Tf;
+ 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_data_ref_put(data));
@@ -346,10 +384,9 @@ main(int argc, char** argv)
scn_args.get_position = box_get_position;
scn_args.nprimitives = box_ntriangles;
scn_args.nvertices = box_nvertices;
- scn_args.trad.temperature = Trad;
- scn_args.trad.reference = Trad;
scn_args.t_range[0] = MMIN(MMIN(Tf, Trad), Tb);
scn_args.t_range[1] = MMAX(MMAX(Tf, Trad), Tb);
+ scn_args.radenv = radenv;
scn_args.context = box_interfaces;
OK(sdis_scene_create(dev, &scn_args, &box_scn));
@@ -359,10 +396,9 @@ main(int argc, char** argv)
scn_args.get_position = square_get_position;
scn_args.nprimitives = square_nsegments;
scn_args.nvertices = square_nvertices;
- scn_args.trad.temperature = Trad;
- scn_args.trad.reference = Trad;
scn_args.t_range[0] = MMIN(MMIN(Tf, Trad), Tb);
scn_args.t_range[1] = MMAX(MMAX(Tf, Trad), Tb);
+ scn_args.radenv = radenv;
scn_args.context = square_interfaces;
OK(sdis_scene_2d_create(dev, &scn_args, &square_scn));
@@ -576,6 +612,7 @@ main(int argc, char** argv)
BA(SOLVE(square_scn, &bound_args, &estimator));
#undef SOLVE
+ OK(sdis_radiative_env_ref_put(radenv));
OK(sdis_scene_ref_put(box_scn));
OK(sdis_scene_ref_put(square_scn));
free_default_device(dev);
@@ -583,4 +620,3 @@ main(int argc, char** argv)
CHK(mem_allocated_size() == 0);
return 0;
}
-
diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c
@@ -253,16 +253,22 @@ interface_get_convection_coef
static double
interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(data != NULL && frag != NULL);
return ((const struct interf*)sdis_data_cget(data))->epsilon;
}
static double
interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(data != NULL && frag != NULL);
return ((const struct interf*)sdis_data_cget(data))->specular_fraction;
}
@@ -284,6 +290,55 @@ interface_get_reference_temperature
}
/*******************************************************************************
+ * Radiative environment
+ ******************************************************************************/
+struct radenv {
+ double temperature; /* [K] */
+ double reference; /* [K] */
+};
+
+static double
+radenv_get_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray;
+ return ((const struct radenv*)sdis_data_cget(data))->temperature;
+}
+
+static double
+radenv_get_reference_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray;
+ return ((const struct radenv*)sdis_data_cget(data))->reference;
+}
+
+static struct sdis_radiative_env*
+create_radenv
+ (struct sdis_device* dev,
+ const double temperature,
+ const double reference)
+{
+ struct sdis_radiative_env_shader shader = SDIS_RADIATIVE_ENV_SHADER_NULL;
+ struct sdis_radiative_env* radenv = NULL;
+ struct sdis_data* data = NULL;
+ struct radenv* env = NULL;
+
+ OK(sdis_data_create(dev, sizeof(struct radenv), ALIGNOF(radenv), NULL, &data));
+ env = sdis_data_get(data);
+ env->temperature = temperature;
+ env->reference = reference;
+
+ shader.temperature = radenv_get_temperature;
+ shader.reference_temperature = radenv_get_reference_temperature;
+ OK(sdis_radiative_env_create(dev, &shader, data, &radenv));
+ OK(sdis_data_ref_put(data));
+ return radenv;
+}
+
+/*******************************************************************************
* Helper functions
******************************************************************************/
static void
@@ -293,7 +348,7 @@ create_solid
struct sdis_medium** solid)
{
struct sdis_data* data = NULL;
- struct solid* solid_param = NULL;
+ struct solid* solid_args = NULL;
struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
CHK(param != NULL);
@@ -302,8 +357,8 @@ create_solid
/* Copy the solid parameters into the Stardis memory space */
OK(sdis_data_create
(dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
- solid_param = sdis_data_get(data);
- memcpy(solid_param, param, sizeof(struct solid));
+ solid_args = sdis_data_get(data);
+ memcpy(solid_args, param, sizeof(struct solid));
/* Setup the solid shader */
solid_shader.calorific_capacity = solid_get_calorific_capacity;
@@ -328,7 +383,7 @@ create_fluid
struct sdis_medium** fluid)
{
struct sdis_data* data = NULL;
- struct fluid* fluid_param = NULL;
+ struct fluid* fluid_args = NULL;
struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
CHK(param != NULL);
@@ -337,8 +392,8 @@ create_fluid
/* Copy the fluid parameters into the Stardis memory space */
OK(sdis_data_create
(dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
- fluid_param = sdis_data_get(data);
- memcpy(fluid_param, param, sizeof(struct fluid));
+ fluid_args = sdis_data_get(data);
+ memcpy(fluid_args, param, sizeof(struct fluid));
/* Setup the fluid shader */
fluid_shader.calorific_capacity = fluid_get_calorific_capacity;
@@ -362,7 +417,7 @@ create_interface
struct sdis_interface** interf)
{
struct sdis_data* data = NULL;
- struct interf* interface_param = NULL;
+ struct interf* interface_args = NULL;
struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
CHK(mdm_front != NULL);
@@ -373,8 +428,8 @@ create_interface
/* Copy the interface parameters into the Stardis memory space */
OK(sdis_data_create
(dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data));
- interface_param = sdis_data_get(data);
- memcpy(interface_param, param, sizeof(struct interf));
+ interface_args = sdis_data_get(data);
+ memcpy(interface_args, param, sizeof(struct interf));
/* Setup the interface shader */
interface_shader.convection_coef = interface_get_convection_coef;
@@ -383,7 +438,7 @@ create_interface
if(sdis_medium_get_type(mdm_front) == SDIS_FLUID) {
interface_shader.front.emissivity = interface_get_emissivity;
interface_shader.front.specular_fraction = interface_get_specular_fraction;
- interface_shader.front.reference_temperature =
+ interface_shader.front.reference_temperature =
interface_get_reference_temperature;
}
if(sdis_medium_get_type(mdm_back) == SDIS_FLUID) {
@@ -555,17 +610,17 @@ main(int argc, char** argv)
struct sdis_medium* fluid1 = NULL;
struct sdis_interface* interf0 = NULL;
struct sdis_interface* interf1 = NULL;
+ struct sdis_radiative_env* radenv = NULL;
struct sdis_scene* scn = NULL;
struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
struct sdis_solve_camera_args solve_args = SDIS_SOLVE_CAMERA_ARGS_DEFAULT;
- struct sdis_ambient_radiative_temperature trad =
- SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
struct ssp_rng* rng = NULL;
struct ssp_rng* rng_state = NULL;
- struct fluid fluid_param = FLUID_NULL;
- struct solid solid_param = SOLID_NULL;
- struct interf interface_param = INTERF_NULL;
- struct fluid* pfluid_param = NULL;
+ struct fluid fluid_args = FLUID_NULL;
+ struct solid solid_args = SOLID_NULL;
+ struct interf interface_args = INTERF_NULL;
+ struct fluid* pfluid_args = NULL;
+ struct radenv* pradenv_args = NULL;
size_t ntris, npos;
size_t nreals, nfails;
size_t definition[2];
@@ -577,40 +632,43 @@ main(int argc, char** argv)
create_default_device(&argc, &argv, &is_master_process, &dev);
+ radenv = create_radenv(dev, 300, 300);
+ pradenv_args = sdis_data_get(sdis_radiative_env_get_data(radenv));
+
/* Create the fluid0 */
- fluid_param.temperature = 350;
- fluid_param.rho = 0;
- fluid_param.cp = 0;
- create_fluid(dev, &fluid_param, &fluid0);
+ fluid_args.temperature = 350;
+ fluid_args.rho = 0;
+ fluid_args.cp = 0;
+ create_fluid(dev, &fluid_args, &fluid0);
/* Create the fluid1 */
- fluid_param.temperature = 300;
- fluid_param.rho = 0;
- fluid_param.cp = 0;
- create_fluid(dev, &fluid_param, &fluid1);
+ fluid_args.temperature = 300;
+ fluid_args.rho = 0;
+ fluid_args.cp = 0;
+ create_fluid(dev, &fluid_args, &fluid1);
/* Create the solid medium */
- solid_param.cp = 1.0;
- solid_param.lambda = 0.1;
- solid_param.rho = 1.0;
- solid_param.delta = 1.0/20.0;
- solid_param.temperature = SDIS_TEMPERATURE_NONE;
- create_solid(dev, &solid_param, &solid);
+ solid_args.cp = 1.0;
+ solid_args.lambda = 0.1;
+ solid_args.rho = 1.0;
+ solid_args.delta = 1.0/20.0;
+ solid_args.temperature = SDIS_TEMPERATURE_NONE;
+ create_solid(dev, &solid_args, &solid);
/* Create the fluid0/solid interface */
- interface_param.hc = 1;
- interface_param.epsilon = 0;
- interface_param.specular_fraction = 0;
- interface_param.temperature = SDIS_TEMPERATURE_NONE;
- create_interface(dev, solid, fluid0, &interface_param, &interf0);
+ interface_args.hc = 1;
+ interface_args.epsilon = 0;
+ interface_args.specular_fraction = 0;
+ interface_args.temperature = SDIS_TEMPERATURE_NONE;
+ create_interface(dev, solid, fluid0, &interface_args, &interf0);
/* Create the fluid1/solid interface */
- interface_param.hc = 0.1;
- interface_param.epsilon = 1;
- interface_param.specular_fraction = 1;
- interface_param.temperature = SDIS_TEMPERATURE_NONE;
- interface_param.reference_temperature = 300;
- create_interface(dev, fluid1, solid, &interface_param, &interf1);
+ interface_args.hc = 0.1;
+ interface_args.epsilon = 1;
+ interface_args.specular_fraction = 1;
+ interface_args.temperature = SDIS_TEMPERATURE_NONE;
+ interface_args.reference_temperature = 300;
+ create_interface(dev, fluid1, solid, &interface_args, &interf1);
/* Setup the cube geometry */
OK(s3dut_create_cuboid(NULL, 2, 2, 2, &msh));
@@ -634,10 +692,9 @@ main(int argc, char** argv)
scn_args.get_position = geometry_get_position;
scn_args.nprimitives = ntris;
scn_args.nvertices = npos;
- scn_args.trad.temperature = 300;
- scn_args.trad.reference = 300;
scn_args.t_range[0] = 300;
scn_args.t_range[1] = 350;
+ scn_args.radenv = radenv;
scn_args.context = &geom;
OK(sdis_scene_create(dev, &scn_args, &scn));
@@ -672,12 +729,9 @@ main(int argc, char** argv)
solve_args.cam = NULL;
BA(sdis_solve_camera(scn, &solve_args, &buf));
solve_args.cam = cam;
- OK(sdis_scene_get_ambient_radiative_temperature(scn, &trad));
- trad.temperature = SDIS_TEMPERATURE_NONE;
- OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
+ pradenv_args->temperature = SDIS_TEMPERATURE_NONE;
BA(sdis_solve_camera(scn, &solve_args, &buf));
- trad.temperature = 300;
- OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
+ pradenv_args->temperature = 300;
solve_args.time_range[0] = solve_args.time_range[1] = -1;
BA(sdis_solve_camera(scn, &solve_args, &buf));
solve_args.time_range[0] = 1;
@@ -765,8 +819,8 @@ main(int argc, char** argv)
solve_args.rng_state = SDIS_SOLVE_CAMERA_ARGS_DEFAULT.rng_state;
solve_args.rng_type = SDIS_SOLVE_CAMERA_ARGS_DEFAULT.rng_type;
- pfluid_param = sdis_data_get(sdis_medium_get_data(fluid1));
- pfluid_param->temperature = SDIS_TEMPERATURE_NONE;
+ pfluid_args = sdis_data_get(sdis_medium_get_data(fluid1));
+ pfluid_args->temperature = SDIS_TEMPERATURE_NONE;
/* Check simulation error handling */
BA(sdis_solve_camera(scn, &solve_args, &buf));
@@ -777,6 +831,7 @@ main(int argc, char** argv)
OK(sdis_medium_ref_put(solid));
OK(sdis_medium_ref_put(fluid0));
OK(sdis_medium_ref_put(fluid1));
+ OK(sdis_radiative_env_ref_put(radenv));
OK(sdis_scene_ref_put(scn));
OK(sdis_camera_ref_put(cam));
OK(sdis_interface_ref_put(interf0));
diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c
@@ -178,16 +178,22 @@ interface_get_convection_coef
static double
interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(data != NULL && frag != NULL);
return ((const struct interf*)sdis_data_cget(data))->epsilon;
}
static double
interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(data != NULL && frag != NULL);
return ((const struct interf*)sdis_data_cget(data))->specular_fraction;
}
diff --git a/src/test_sdis_solve_medium_2d.c b/src/test_sdis_solve_medium_2d.c
@@ -166,16 +166,22 @@ interface_get_convection_coef
static double
interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(data != NULL && frag != NULL);
return ((const struct interf*)sdis_data_cget(data))->epsilon;
}
static double
interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(data != NULL && frag != NULL);
return ((const struct interf*)sdis_data_cget(data))->specular_fraction;
}
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -156,16 +156,22 @@ interface_get_convection_coef
static double
interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(data != NULL && frag != NULL);
return ((const struct interf*)sdis_data_cget(data))->epsilon;
}
static double
interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
{
+ (void)source_id;
CHK(data != NULL && frag != NULL);
return ((const struct interf*)sdis_data_cget(data))->specular_fraction;
}
@@ -179,6 +185,55 @@ interface_get_reference_temperature
}
/*******************************************************************************
+ * Radiative environment
+ ******************************************************************************/
+struct radenv {
+ double temperature; /* [K] */
+ double reference; /* [K] */
+};
+
+static double
+radenv_get_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray;
+ return ((const struct radenv*)sdis_data_cget(data))->temperature;
+}
+
+static double
+radenv_get_reference_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray;
+ return ((const struct radenv*)sdis_data_cget(data))->reference;
+}
+
+static struct sdis_radiative_env*
+create_radenv
+ (struct sdis_device* dev,
+ const double temperature, /* [K] */
+ const double reference) /* [K] */
+{
+ struct sdis_radiative_env_shader shader = SDIS_RADIATIVE_ENV_SHADER_NULL;
+ struct sdis_radiative_env* radenv = NULL;
+ struct sdis_data* data = NULL;
+ struct radenv* radenv_args = NULL;
+
+ OK(sdis_data_create(dev, sizeof(struct radenv), ALIGNOF(radenv), NULL, &data));
+ radenv_args = sdis_data_get(data);
+ radenv_args->temperature = temperature;
+ radenv_args->reference = reference;
+
+ shader.temperature = radenv_get_temperature;
+ shader.reference_temperature = radenv_get_reference_temperature;
+ OK(sdis_radiative_env_create(dev, &shader, data, &radenv));
+ OK(sdis_data_ref_put(data));
+ return radenv;
+}
+
+/*******************************************************************************
* Helper functions
******************************************************************************/
struct dump_path_context {
@@ -274,6 +329,7 @@ main(int argc, char** argv)
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid = NULL;
struct sdis_interface* interf = NULL;
+ struct sdis_radiative_env* radenv = NULL;
struct sdis_scene* scn = NULL;
struct sdis_data* data = NULL;
struct sdis_estimator* estimator = NULL;
@@ -289,13 +345,12 @@ main(int argc, char** argv)
struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
- struct sdis_ambient_radiative_temperature trad =
- SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
struct dump_path_context dump_ctx = DUMP_PATH_CONTEXT_NULL;
struct context ctx;
- struct fluid* fluid_param;
- struct solid* solid_param;
- struct interf* interface_param;
+ struct radenv* radenv_args;
+ struct fluid* fluid_args;
+ struct solid* solid_args;
+ struct interf* interface_args;
struct ssp_rng* rng_state = NULL;
enum sdis_estimator_type type;
FILE* stream = NULL;
@@ -315,8 +370,8 @@ main(int argc, char** argv)
/* Create the fluid medium */
OK(sdis_data_create
(dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
- fluid_param = sdis_data_get(data);
- fluid_param->temperature = 300;
+ fluid_args = sdis_data_get(data);
+ fluid_args->temperature = 300;
fluid_shader.temperature = fluid_get_temperature;
OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid));
OK(sdis_data_ref_put(data));
@@ -324,12 +379,12 @@ main(int argc, char** argv)
/* Create the solid medium */
OK(sdis_data_create
(dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
- solid_param = sdis_data_get(data);
- solid_param->cp = 1.0;
- solid_param->lambda = 0.1;
- solid_param->rho = 1.0;
- solid_param->delta = 1.0/20.0;
- solid_param->temperature = SDIS_TEMPERATURE_NONE; /* Unknown temperature */
+ solid_args = sdis_data_get(data);
+ solid_args->cp = 1.0;
+ solid_args->lambda = 0.1;
+ solid_args->rho = 1.0;
+ solid_args->delta = 1.0/20.0;
+ solid_args->temperature = SDIS_TEMPERATURE_NONE; /* Unknown temperature */
solid_shader.calorific_capacity = solid_get_calorific_capacity;
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
@@ -341,10 +396,10 @@ main(int argc, char** argv)
/* Create the solid/fluid interface */
OK(sdis_data_create(dev, sizeof(struct interf),
ALIGNOF(struct interf), NULL, &data));
- interface_param = sdis_data_get(data);
- interface_param->hc = 0.5;
- interface_param->epsilon = 0;
- interface_param->specular_fraction = 0;
+ interface_args = sdis_data_get(data);
+ interface_args->hc = 0.5;
+ interface_args->epsilon = 0;
+ interface_args->specular_fraction = 0;
interface_shader.convection_coef = interface_get_convection_coef;
interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL;
interface_shader.back.temperature = NULL;
@@ -359,6 +414,10 @@ main(int argc, char** argv)
OK(sdis_medium_ref_put(solid));
OK(sdis_medium_ref_put(fluid));
+ /* Create the radiative environment */
+ radenv = create_radenv(dev, SDIS_TEMPERATURE_NONE, SDIS_TEMPERATURE_NONE);
+ radenv_args = sdis_data_get(sdis_radiative_env_get_data(radenv));
+
/* Create the scene */
ctx.positions = box_vertices;
ctx.indices = box_indices;
@@ -368,6 +427,7 @@ main(int argc, char** argv)
scn_args.get_position = get_position;
scn_args.nprimitives = box_ntriangles;
scn_args.nvertices = box_nvertices;
+ scn_args.radenv = radenv;
scn_args.context = &ctx;
OK(sdis_scene_create(dev, &scn_args, &scn));
@@ -442,7 +502,7 @@ main(int argc, char** argv)
ref = 300;
printf("Temperature at (%g, %g, %g) with Tfluid=%g = %g ~ %g +/- %g\n",
- SPLIT3(solve_args.position), fluid_param->temperature, ref, T.E, T.SE);
+ SPLIT3(solve_args.position), fluid_args->temperature, ref, T.E, T.SE);
printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
@@ -457,10 +517,10 @@ main(int argc, char** argv)
OK(sdis_estimator_ref_put(estimator));
/* The external fluid cannot have an unknown temperature */
- fluid_param->temperature = SDIS_TEMPERATURE_NONE;
+ fluid_args->temperature = SDIS_TEMPERATURE_NONE;
BA(sdis_solve_probe(scn, &solve_args, &estimator));
- fluid_param->temperature = 300;
+ fluid_args->temperature = 300;
OK(sdis_solve_probe(scn, &solve_args, &estimator));
BA(sdis_solve_probe_green_function(NULL, &solve_args, &green));
@@ -484,7 +544,7 @@ main(int argc, char** argv)
printf("\n");
/* Check same green used at a different temperature */
- fluid_param->temperature = 500;
+ fluid_args->temperature = 500;
OK(sdis_solve_probe(scn, &solve_args, &estimator));
OK(sdis_estimator_get_realisation_count(estimator, &nreals));
@@ -494,7 +554,7 @@ main(int argc, char** argv)
ref = 500;
printf("Temperature at (%g, %g, %g) with Tfluid=%g = %g ~ %g +/- %g\n",
- SPLIT3(solve_args.position), fluid_param->temperature, ref, T.E, T.SE);
+ SPLIT3(solve_args.position), fluid_args->temperature, ref, T.E, T.SE);
printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
@@ -571,10 +631,10 @@ main(int argc, char** argv)
solve_args.register_paths = SDIS_HEAT_PATH_ALL;
/* Check simulation error handling when paths are registered */
- fluid_param->temperature = SDIS_TEMPERATURE_NONE;
+ fluid_args->temperature = SDIS_TEMPERATURE_NONE;
BA(sdis_solve_probe(scn, &solve_args, &estimator));
- fluid_param->temperature = 300;
+ fluid_args->temperature = 300;
OK(sdis_solve_probe(scn, &solve_args, &estimator));
OK(sdis_estimator_get_paths_count(estimator, &n));
CHK(n == N_dump);
@@ -595,14 +655,14 @@ main(int argc, char** argv)
/* Green and ambient radiative temperature */
solve_args.nrealisations = N;
- trad.temperature = trad.reference = 300;
+ radenv_args->temperature = 300;
+ radenv_args->reference = 300;
t_range[0] = 300;
t_range[1] = 300;
- OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
OK(sdis_scene_set_temperature_range(scn, t_range));
- interface_param->epsilon = 1;
- interface_param->reference_temperature = 300;
+ interface_args->epsilon = 1;
+ interface_args->reference_temperature = 300;
OK(sdis_solve_probe(scn, &solve_args, &estimator));
OK(sdis_solve_probe_green_function(scn, &solve_args, &green));
@@ -615,10 +675,9 @@ main(int argc, char** argv)
OK(sdis_estimator_ref_put(estimator2));
/* Check same green used at different ambient radiative temperature */
- trad.temperature = 600;
+ radenv_args->temperature = 300;
t_range[0] = 300;
t_range[1] = 600;
- OK(sdis_scene_set_ambient_radiative_temperature(scn, &trad));
OK(sdis_scene_set_temperature_range(scn, t_range));
OK(sdis_solve_probe(scn, &solve_args, &estimator));
@@ -630,6 +689,7 @@ main(int argc, char** argv)
OK(sdis_estimator_ref_put(estimator));
OK(sdis_estimator_ref_put(estimator2));
OK(sdis_green_function_ref_put(green));
+ OK(sdis_radiative_env_ref_put(radenv));
OK(sdis_scene_ref_put(scn));
OK(sdis_device_ref_put(dev));
diff --git a/src/test_sdis_source.c b/src/test_sdis_source.c
@@ -44,6 +44,7 @@ check_spherical_source(struct sdis_device* dev)
struct sdis_spherical_source_create_args args =
SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL;
struct sdis_source* src = NULL;
+ struct sdis_source* src2 = NULL;
struct sdis_data* data = NULL;
/* Create a data to check its memory management */
@@ -71,7 +72,10 @@ check_spherical_source(struct sdis_device* dev)
args.data = NULL;
OK(sdis_spherical_source_create(dev, &args, &src));
+ OK(sdis_spherical_source_create(dev, &args, &src2));
+ CHK(sdis_source_get_id(src) != sdis_source_get_id(src2));
OK(sdis_source_ref_put(src));
+ OK(sdis_source_ref_put(src2));
args.position = NULL;
BA(sdis_spherical_source_create(dev, &args, &src));
diff --git a/src/test_sdis_unstationary_atm.c b/src/test_sdis_unstationary_atm.c
@@ -1,918 +0,0 @@
-/* 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 "sdis.h"
-#include "test_sdis_utils.h"
-
-#include <rsys/clock_time.h>
-#include <rsys/mem_allocator.h>
-#include <rsys/double3.h>
-#include <rsys/math.h>
-#include <star/ssp.h>
-
-/*
- * The physical configuration is the following: a slab of fluid with known
- * thermophysical properties but unknown temperature is located between a
- * "ground" and a slab of solid, with also a unknown temperature profile. On
- * the other side of the solid slab, is an "atmosphere" with known temperature,
- * and known radiative temperature.
- *
- * Solving the system means: finding the temperature of the ground, of the
- * fluid, of the boundaries, and also the temperature inside the solid, at
- * various locations (the 1D slab is discretized in order to obtain the
- * reference)
- *
- * The reference for this system comes from a numerical method and is not
- * analytic. Thus the compliance test MC VS reference is not the usual |MC -
- * ref| <= 3*sigma but is |MC -ref| <= (Tmax -Tmin) * 0.01.
- *
- * 3D 2D
- *
- * /////////////// ///////////////
- * +-----------+-+ +-----------+-+
- * /' / /| | | |
- * +-----------+-+ | HA _\ <---- TR | _\ | | HA _\ <---- TR
- * | | _\ | | | / / <---- TR TG| HG / / HC | | / / <---- TR
- * | |HG / / HC| | | TA \__/ <---- TR | \__/ | | TA \__/ <---- TR
- * TG| | \__/ | | | | | |
- * | +.........|.|.+ +-----------+-+
- * |, |/|/ /////H///////E/
- * +-----------+-+
- * //////H//////E/
- */
-
-#define XH 3
-#define XHpE 3.2
-#define XE (XHpE - XH)
-
-#define T0_SOLID 300
-#define T0_FLUID 300
-
-#define N 10000 /* #realisations */
-
-#define TG 310
-#define HG 400
-
-#define HC 400
-
-#define TA 290
-#define HA 400
-#define TR 260
-
-#define TMAX (MMAX(MMAX(MMAX(T0_FLUID, T0_SOLID), MMAX(TG, TA)), TR))
-#define TMIN (MMIN(MMIN(MMIN(T0_FLUID, T0_SOLID), MMIN(TG, TA)), TR))
-#define EPS ((TMAX-TMIN)*0.01)
-
-/* hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon
- * Tref = (hr / (4 * 5.6696e-8 * epsilon)) ^ 1/3, hr = 6 */
-#define TREF 297.974852286
-
-#define RHO_F 1.3
-#define CP_F 1005
-#define RHO_S 2400
-#define CP_S 800
-#define LAMBDA 0.6
-
-#define X_PROBE (XH + 0.2 * XE)
-
-#define DELTA (XE/40.0)
-
-/*******************************************************************************
- * Box geometry
- ******************************************************************************/
-static const double model3d_vertices[12/*#vertices*/*3/*#coords per vertex*/] = {
- 0, 0, 0,
- XH, 0, 0,
- XHpE, 0, 0,
- 0, XHpE, 0,
- XH, XHpE, 0,
- XHpE, XHpE, 0,
- 0, 0, XHpE,
- XH, 0, XHpE,
- XHpE, 0, XHpE,
- 0, XHpE, XHpE,
- XH, XHpE, XHpE,
- XHpE, XHpE, XHpE
-};
-static const size_t model3d_nvertices = sizeof(model3d_vertices)/(sizeof(double)*3);
-
-/* The following array lists the indices toward the 3D vertices of each
- * triangle.
- * ,3---,4---,5 ,3----4----5 ,4
- * ,' | ,' | ,'/| ,'/| \ | \ | ,'/|
- * 9----10---11 / | 9' / | \ | \ | 10 / | Y
- * |', |', | / ,2 | / ,0---,1---,2 | / ,1 |
- * | ',| ',|/,' |/,' | ,' | ,' |/,' o--X
- * 6----7----8' 6----7'---8' 7 /
- * Front, right Back, left and Internal Z
- * and Top faces bottom faces face */
-static const size_t model3d_indices[22/*#triangles*/*3/*#indices per triangle*/] = {
- 0, 3, 1, 1, 3, 4, 1, 4, 2, 2, 4, 5, /* -Z */
- 0, 6, 3, 3, 6, 9, /* -X */
- 6, 7, 9, 9, 7, 10, 7, 8, 10, 10, 8, 11, /* +Z */
- 5, 11, 8, 8, 2, 5, /* +X */
- 3, 9, 10, 10, 4, 3, 4, 10, 11, 11, 5, 4, /* +Y */
- 0, 1, 7, 7, 6, 0, 1, 2, 8, 8, 7, 1, /* -Y */
- 4, 10, 7, 7, 1, 4 /* Inside */
-};
-static const size_t model3d_ntriangles = sizeof(model3d_indices)/(sizeof(size_t)*3);
-
-static INLINE void
-model3d_get_indices(const size_t itri, size_t ids[3], void* context)
-{
- (void)context;
- CHK(ids);
- CHK(itri < model3d_ntriangles);
- ids[0] = model3d_indices[itri * 3 + 0];
- ids[1] = model3d_indices[itri * 3 + 1];
- ids[2] = model3d_indices[itri * 3 + 2];
-}
-
-static INLINE void
-model3d_get_position(const size_t ivert, double pos[3], void* context)
-{
- (void)context;
- CHK(pos);
- CHK(ivert < model3d_nvertices);
- pos[0] = model3d_vertices[ivert * 3 + 0];
- pos[1] = model3d_vertices[ivert * 3 + 1];
- pos[2] = model3d_vertices[ivert * 3 + 2];
-}
-
-static INLINE void
-model3d_get_interface(const size_t itri, struct sdis_interface** bound, void* context)
-{
- struct sdis_interface** interfaces = context;
- CHK(context && bound);
- CHK(itri < model3d_ntriangles);
- *bound = interfaces[itri];
-}
-
-/*******************************************************************************
- * Square geometry
- ******************************************************************************/
-static const double model2d_vertices[6/*#vertices*/ * 2/*#coords per vertex*/] = {
- XHpE, 0,
- XH, 0,
- 0, 0,
- 0, XHpE,
- XH, XHpE,
- XHpE, XHpE
-};
-static const size_t model2d_nvertices = sizeof(model2d_vertices)/(sizeof(double)*2);
-
-static const size_t model2d_indices[7/*#segments*/ * 2/*#indices per segment*/] = {
- 0, 1, 1, 2, /* Bottom */
- 2, 3, /* Left */
- 3, 4, 4, 5, /* Top */
- 5, 0, /* Right */
- 4, 1 /* Inside */
-};
-static const size_t model2d_nsegments = sizeof(model2d_indices)/(sizeof(size_t)*2);
-
-static INLINE void
-model2d_get_indices(const size_t iseg, size_t ids[2], void* context)
-{
- (void)context;
- CHK(ids);
- CHK(iseg < model2d_nsegments);
- ids[0] = model2d_indices[iseg * 2 + 0];
- ids[1] = model2d_indices[iseg * 2 + 1];
-}
-
-static INLINE void
-model2d_get_position(const size_t ivert, double pos[2], void* context)
-{
- (void)context;
- CHK(pos);
- CHK(ivert < model2d_nvertices);
- pos[0] = model2d_vertices[ivert * 2 + 0];
- pos[1] = model2d_vertices[ivert * 2 + 1];
-}
-
-static INLINE void
-model2d_get_interface
-(const size_t iseg, struct sdis_interface** bound, void* context)
-{
- struct sdis_interface** interfaces = context;
- CHK(context && bound);
- CHK(iseg < model2d_nsegments);
- *bound = interfaces[iseg];
-}
-
-/*******************************************************************************
- * Media
- ******************************************************************************/
-struct solid {
- double lambda;
- double rho;
- double cp;
- double delta;
- double temperature;
- double t0;
-};
-
-static double
-solid_get_calorific_capacity
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- struct solid* solid;
- CHK(vtx && data);
- solid = ((struct solid*)sdis_data_cget(data));
- return solid->cp;
-}
-
-static double
-solid_get_thermal_conductivity
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- struct solid* solid;
- CHK(vtx && data);
- solid = ((struct solid*)sdis_data_cget(data));
- return solid->lambda;
-}
-
-static double
-solid_get_volumic_mass
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- struct solid* solid;
- CHK(vtx && data);
- solid = ((struct solid*)sdis_data_cget(data));
- return solid->rho;
-}
-
-static double
-solid_get_delta
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- struct solid* solid;
- CHK(vtx && data);
- solid = ((struct solid*)sdis_data_cget(data));
- return solid->delta;
-}
-
-static double
-solid_get_temperature
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- struct solid* solid;
- CHK(vtx && data);
- solid = ((struct solid*)sdis_data_cget(data));
- if(vtx->time <= solid->t0)
- return solid->temperature;
- return SDIS_TEMPERATURE_NONE;
-}
-
-struct fluid {
- double rho;
- double cp;
- double t0;
- double temperature;
-};
-
-static double
-fluid_get_temperature
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- struct fluid* fluid;
- CHK(vtx && data);
- fluid = ((struct fluid*)sdis_data_cget(data));
- if(vtx->time <= fluid->t0)
- return fluid->temperature;
- return SDIS_TEMPERATURE_NONE;
-}
-
-static double
-fluid_get_volumic_mass
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- struct fluid* fluid;
- CHK(vtx && data);
- fluid = ((struct fluid*)sdis_data_cget(data));
- return fluid->rho;
-}
-
-static double
-fluid_get_calorific_capacity
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- struct fluid* fluid;
- CHK(vtx && data);
- fluid = ((struct fluid*)sdis_data_cget(data));
- return fluid->cp;
-}
-
-/*******************************************************************************
- * Interfaces
- ******************************************************************************/
-struct interf {
- double temperature;
- double emissivity;
- double h;
- double Tref;
-};
-
-static double
-interface_get_temperature
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- const struct interf* interf;
- CHK(frag && data);
- interf = sdis_data_cget(data);
- return interf->temperature;
-}
-
-static double
-interface_get_convection_coef
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- const struct interf* interf;
- CHK(frag && data);
- interf = sdis_data_cget(data);
- return interf->h;
-}
-
-static double
-interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- const struct interf* interf;
- CHK(frag && data);
- interf = sdis_data_cget(data);
- return interf->emissivity;
-}
-
-static double
-interface_get_Tref
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- const struct interf* interf;
- CHK(frag && data);
- interf = sdis_data_cget(data);
- return interf->Tref;
-}
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static void
-create_interface
- (struct sdis_device* dev,
- struct sdis_medium* front,
- struct sdis_medium* back,
- const struct interf* interf,
- struct sdis_interface** out_interf)
-{
- struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
- struct sdis_data* data = NULL;
-
- CHK(interf != NULL);
-
- shader.front.temperature = interface_get_temperature;
- shader.back.temperature = interface_get_temperature;
- if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) {
- shader.convection_coef = interface_get_convection_coef;
- shader.convection_coef_upper_bound = interf->h;
- }
- if(sdis_medium_get_type(front) == SDIS_FLUID) {
- shader.front.emissivity = interface_get_emissivity;
- shader.front.reference_temperature = interface_get_Tref;
- }
- if(sdis_medium_get_type(back) == SDIS_FLUID) {
- shader.back.emissivity = interface_get_emissivity;
- shader.back.reference_temperature = interface_get_Tref;
- }
-
- OK(sdis_data_create(dev, sizeof(struct interf), ALIGNOF(struct interf),
- NULL, &data));
- *((struct interf*)sdis_data_get(data)) = *interf;
-
- OK(sdis_interface_create(dev, front, back, &shader, data, out_interf));
- OK(sdis_data_ref_put(data));
-}
-
-static void
-solve_tbound1
- (struct sdis_scene* scn,
- struct ssp_rng* rng)
-{
- char dump[128];
- struct time t0, t1;
- struct sdis_estimator* estimator;
- struct sdis_solve_probe_boundary_args solve_args
- = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT;
- struct sdis_mc T = SDIS_MC_NULL;
- struct sdis_mc time = SDIS_MC_NULL;
- size_t nreals;
- size_t nfails;
- enum sdis_scene_dimension dim;
- const double t[] = { 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
- const double ref[sizeof(t) / sizeof(*t)] = {
- 290.046375, 289.903935, 289.840490, 289.802690, 289.777215, 289.759034,
- 289.745710, 289.735826, 289.728448, 289.722921
- };
- const int nsimuls = sizeof(t) / sizeof(*t);
- int isimul;
- ASSERT(scn && rng);
-
- OK(sdis_scene_get_dimension(scn, &dim));
-
- solve_args.nrealisations = N;
- solve_args.side = SDIS_FRONT;
- FOR_EACH(isimul, 0, nsimuls) {
- solve_args.time_range[0] = solve_args.time_range[1] = t[isimul];
- if(dim == SDIS_SCENE_2D) {
- solve_args.iprim = model2d_nsegments - 1;
- solve_args.uv[0] = ssp_rng_uniform_double(rng, 0, 1);
- } else {
- double u = ssp_rng_uniform_double(rng, 0, 1);
- double v = ssp_rng_uniform_double(rng, 0, 1);
- double w = ssp_rng_uniform_double(rng, 0, 1);
- double x = 1 / (u + v + w);
- solve_args.iprim = (isimul % 2) ? 10 : 11; /* +X face */
- solve_args.uv[0] = u * x;
- solve_args.uv[1] = v * x;
- }
-
- time_current(&t0);
- OK(sdis_solve_probe_boundary(scn, &solve_args, &estimator));
- time_sub(&t0, time_current(&t1), &t0);
- time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
-
- OK(sdis_estimator_get_realisation_count(estimator, &nreals));
- OK(sdis_estimator_get_failure_count(estimator, &nfails));
- OK(sdis_estimator_get_temperature(estimator, &T));
- OK(sdis_estimator_get_realisation_time(estimator, &time));
-
- switch(dim) {
- case SDIS_SCENE_2D:
- printf("Unstationary temperature at (%lu/%g) at t=%g = %g ~ %g +/- %g\n",
- (unsigned long)solve_args.iprim, solve_args.uv[0], t[isimul],
- ref[isimul], T.E, T.SE);
- break;
- case SDIS_SCENE_3D:
- printf("Unstationary temperature at (%lu/%g,%g) at t=%g = %g ~ %g +/- %g\n",
- (unsigned long)solve_args.iprim, SPLIT2(solve_args.uv), t[isimul],
- ref[isimul], T.E, T.SE);
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- printf("Elapsed time = %s\n", dump);
- printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
-
- CHK(eq_eps(T.E, ref[isimul], EPS));
- /*CHK(eq_eps(T.E, ref[isimul], T.SE*3));*/
-
- OK(sdis_estimator_ref_put(estimator));
- }
-}
-
-static void
-solve_tbound2
- (struct sdis_scene* scn,
- struct ssp_rng* rng)
-{
- char dump[128];
- struct time t0, t1;
- struct sdis_estimator* estimator;
- struct sdis_solve_probe_boundary_args solve_args
- = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT;
- struct sdis_mc T = SDIS_MC_NULL;
- struct sdis_mc time = SDIS_MC_NULL;
- size_t nreals;
- size_t nfails;
- enum sdis_scene_dimension dim;
- const double t[] = { 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
- const double ref[sizeof(t) / sizeof(*t)] = {
- 309.08032, 309.34626, 309.46525, 309.53625, 309.58408, 309.618121,
- 309.642928, 309.661167, 309.674614, 309.684524
- };
- const int nsimuls = sizeof(t) / sizeof(*t);
- int isimul;
- ASSERT(scn && rng);
-
- OK(sdis_scene_get_dimension(scn, &dim));
-
- solve_args.nrealisations = N;
- solve_args.side = SDIS_FRONT;
- FOR_EACH(isimul, 0, nsimuls) {
- solve_args.time_range[0] = solve_args.time_range[1] = t[isimul];
- if(dim == SDIS_SCENE_2D) {
- solve_args.iprim = model2d_nsegments - 1;
- solve_args.uv[0] = ssp_rng_uniform_double(rng, 0, 1);
- } else {
- double u = ssp_rng_uniform_double(rng, 0, 1);
- double v = ssp_rng_uniform_double(rng, 0, 1);
- double w = ssp_rng_uniform_double(rng, 0, 1);
- double x = 1 / (u + v + w);
- solve_args.iprim = (isimul % 2) ? 20 : 21; /* Internal face */
- solve_args.uv[0] = u * x;
- solve_args.uv[1] = v * x;
- }
-
- time_current(&t0);
- OK(sdis_solve_probe_boundary(scn, &solve_args, &estimator));
- time_sub(&t0, time_current(&t1), &t0);
- time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
-
- OK(sdis_estimator_get_realisation_count(estimator, &nreals));
- OK(sdis_estimator_get_failure_count(estimator, &nfails));
- OK(sdis_estimator_get_temperature(estimator, &T));
- OK(sdis_estimator_get_realisation_time(estimator, &time));
-
- switch(dim) {
- case SDIS_SCENE_2D:
- printf("Unstationary temperature at (%lu/%g) at t=%g = %g ~ %g +/- %g\n",
- (unsigned long)solve_args.iprim, solve_args.uv[0], t[isimul],
- ref[isimul], T.E, T.SE);
- break;
- case SDIS_SCENE_3D:
- printf("Unstationary temperature at (%lu/%g,%g) at t=%g = %g ~ %g +/- %g\n",
- (unsigned long)solve_args.iprim, SPLIT2(solve_args.uv), t[isimul],
- ref[isimul], T.E, T.SE);
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- printf("Elapsed time = %s\n", dump);
- printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
-
- CHK(nfails + nreals == N);
- CHK(nfails <= N/1000);
- CHK(eq_eps(T.E, ref[isimul], EPS));
- /*CHK(eq_eps(T.E, ref[isimul], T.SE*3));*/
-
- OK(sdis_estimator_ref_put(estimator));
- }
-}
-
-static void
-solve_tsolid
- (struct sdis_scene* scn,
- struct ssp_rng* rng)
-{
- char dump[128];
- struct time t0, t1;
- struct sdis_estimator* estimator;
- struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
- struct sdis_mc T = SDIS_MC_NULL;
- struct sdis_mc time = SDIS_MC_NULL;
- size_t nreals;
- size_t nfails;
- enum sdis_scene_dimension dim;
- const double t[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
- const double ref[sizeof(t) / sizeof(*t)] = {
- 300, 300.87408, 302.25832, 303.22164, 303.89954, 304.39030, 304.75041,
- 305.01595, 305.21193, 305.35641, 305.46271
- };
- const int nsimuls = sizeof(t) / sizeof(*t);
- int isimul;
- ASSERT(scn && rng);
-
- OK(sdis_scene_get_dimension(scn, &dim));
-
- solve_args.nrealisations = N;
- FOR_EACH(isimul, 0, nsimuls) {
- solve_args.time_range[0] = solve_args.time_range[1] = t[isimul];
- solve_args.position[0] = X_PROBE;
- solve_args.position[1] = ssp_rng_uniform_double(rng, 0.1*XHpE, 0.9*XHpE);
-
- if(dim == SDIS_SCENE_3D)
- solve_args.position[2] = ssp_rng_uniform_double(rng, 0.1*XHpE, 0.9*XHpE);
-
- time_current(&t0);
- OK(sdis_solve_probe(scn, &solve_args, &estimator));
- time_sub(&t0, time_current(&t1), &t0);
- time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
-
- OK(sdis_estimator_get_realisation_count(estimator, &nreals));
- OK(sdis_estimator_get_failure_count(estimator, &nfails));
- OK(sdis_estimator_get_temperature(estimator, &T));
- OK(sdis_estimator_get_realisation_time(estimator, &time));
-
- switch (dim) {
- case SDIS_SCENE_2D:
- printf("Unstationary temperature at (%g,%g) at t=%g = %g ~ %g +/- %g\n",
- SPLIT2(solve_args.position), t[isimul], ref[isimul], T.E, T.SE);
- break;
- case SDIS_SCENE_3D:
- printf("Unstationary temperature at (%g,%g,%g) at t=%g = %g ~ %g +/- %g\n",
- SPLIT3(solve_args.position), t[isimul], ref[isimul], T.E, T.SE);
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- printf("Elapsed time = %s\n", dump);
- printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
-
- CHK(nfails + nreals == N);
- CHK(nfails <= N / 1000);
- CHK(eq_eps(T.E, ref[isimul], EPS));
- /*CHK(eq_eps(T.E, ref[isimul], T.SE*3));*/
-
- OK(sdis_estimator_ref_put(estimator));
- }
-}
-
-static void
-solve_tfluid
- (struct sdis_scene* scn)
-{
- char dump[128];
- struct time t0, t1;
- struct sdis_estimator* estimator;
- struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
- struct sdis_mc T = SDIS_MC_NULL;
- struct sdis_mc time = SDIS_MC_NULL;
- size_t nreals;
- size_t nfails;
- enum sdis_scene_dimension dim;
- double eps;
- const double t[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
- const double ref[sizeof(t) / sizeof(*t)] = {
- 300, 309.53905, 309.67273, 309.73241, 309.76798, 309.79194, 309.80899,
- 309.82141, 309.83055, 309.83728, 309.84224
- };
- const int nsimuls = sizeof(t) / sizeof(*t);
- int isimul;
- ASSERT(scn);
-
- OK(sdis_scene_get_dimension(scn, &dim));
-
- solve_args.nrealisations = N;
- solve_args.position[0] = XH * 0.5;
- solve_args.position[1] = XH * 0.5;
- solve_args.position[2] = XH * 0.5;
- FOR_EACH(isimul, 0, nsimuls) {
- solve_args.time_range[0] = solve_args.time_range[1] = t[isimul];
-
- time_current(&t0);
- OK(sdis_solve_probe(scn, &solve_args, &estimator));
- time_sub(&t0, time_current(&t1), &t0);
- time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
-
- OK(sdis_estimator_get_realisation_count(estimator, &nreals));
- OK(sdis_estimator_get_failure_count(estimator, &nfails));
- OK(sdis_estimator_get_temperature(estimator, &T));
- OK(sdis_estimator_get_realisation_time(estimator, &time));
-
- switch (dim) {
- case SDIS_SCENE_2D:
- printf("Unstationary fluid temperature at t=%g = %g ~ %g +/- %g\n",
- t[isimul], ref[isimul], T.E, T.SE);
- break;
- case SDIS_SCENE_3D:
- printf("Unstationary fluid temperature at t=%g = %g ~ %g +/- %g\n",
- t[isimul], ref[isimul], T.E, T.SE);
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- printf("Elapsed time = %s\n", dump);
- printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
-
- CHK(nfails + nreals == N);
- CHK(nfails <= N / 1000);
-
- eps = EPS;
- CHK(eq_eps(T.E, ref[isimul], eps));
-
- OK(sdis_estimator_ref_put(estimator));
- }
-}
-
-/*******************************************************************************
- * Test
- ******************************************************************************/
-int
-main(int argc, char** argv)
-{
- struct sdis_data* data = NULL;
- struct sdis_device* dev = NULL;
- struct sdis_medium* fluid = NULL;
- struct sdis_medium* fluid_A = NULL;
- struct sdis_medium* solid = NULL;
- struct sdis_medium* dummy_solid = NULL;
- struct sdis_interface* interf_adiabatic_1 = NULL;
- struct sdis_interface* interf_adiabatic_2 = NULL;
- struct sdis_interface* interf_TG = NULL;
- struct sdis_interface* interf_P = NULL;
- struct sdis_interface* interf_TA = NULL;
- struct sdis_scene* box_scn = NULL;
- struct sdis_scene* square_scn = NULL;
- struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
- struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
- struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
- struct sdis_interface* model3d_interfaces[22 /*#triangles*/];
- struct sdis_interface* model2d_interfaces[7/*#segments*/];
- struct interf interf_props;
- struct solid* solid_props = NULL;
- struct fluid* fluid_props = NULL;
- struct ssp_rng* rng = NULL;
- (void)argc, (void)argv;
-
- OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev));
-
- /* Setup the solid shader */
- solid_shader.calorific_capacity = solid_get_calorific_capacity;
- solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
- solid_shader.volumic_mass = solid_get_volumic_mass;
- solid_shader.delta = solid_get_delta;
- solid_shader.temperature = solid_get_temperature;
-
- /* Create the solid media */
- OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data));
- solid_props = sdis_data_get(data);
- solid_props->lambda = LAMBDA;
- solid_props->cp = CP_S;
- solid_props->rho = RHO_S;
- solid_props->delta = DELTA;
- solid_props->t0 = 0;
- solid_props->temperature = T0_SOLID;
- OK(sdis_solid_create(dev, &solid_shader, data, &solid));
- OK(sdis_data_ref_put(data));
-
- /* Create a dummy solid media to be used outside the model */
- OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data));
- solid_props = sdis_data_get(data);
- solid_props->lambda = 0;
- solid_props->cp = 1;
- solid_props->rho = 1;
- solid_props->delta = 1;
- solid_props->t0 = INF;
- solid_props->temperature = SDIS_TEMPERATURE_NONE;
- OK(sdis_solid_create(dev, &solid_shader, data, &dummy_solid));
- OK(sdis_data_ref_put(data));
-
- /* Setup the fluid shader */
- fluid_shader.calorific_capacity = fluid_get_calorific_capacity;
- fluid_shader.volumic_mass = fluid_get_volumic_mass;
- fluid_shader.temperature = fluid_get_temperature;
-
- /* Create the internal fluid media */
- OK(sdis_data_create(dev, sizeof(struct fluid), 16, NULL, &data));
- fluid_props = sdis_data_get(data);
- fluid_props->cp = CP_F;
- fluid_props->rho = RHO_F;
- fluid_props->t0 = 0;
- fluid_props->temperature = T0_FLUID;
- OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid));
- OK(sdis_data_ref_put(data));
-
- /* Create the 'A' fluid media */
- OK(sdis_data_create(dev, sizeof(struct fluid), 16, NULL, &data));
- fluid_props = sdis_data_get(data);
- fluid_props->cp = 1;
- fluid_props->rho = 1;
- fluid_props->t0 = INF;
- fluid_props->temperature = TA;
- OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid_A));
- OK(sdis_data_ref_put(data));
-
- /* Create the adiabatic interfaces */
- interf_props.temperature = SDIS_TEMPERATURE_NONE;
- interf_props.h = 0;
- interf_props.emissivity = 0;
- interf_props.Tref = TREF;
- create_interface(dev, fluid, dummy_solid, &interf_props, &interf_adiabatic_1);
- create_interface(dev, solid, dummy_solid, &interf_props, &interf_adiabatic_2);
-
- /* Create the P interface */
- interf_props.temperature = SDIS_TEMPERATURE_NONE;
- interf_props.h = HC;
- interf_props.emissivity = 1;
- interf_props.Tref = TREF;
- create_interface(dev, fluid, solid, &interf_props, &interf_P);
-
- /* Create the TG interface */
- interf_props.temperature = TG;
- interf_props.h = HG;
- interf_props.emissivity = 1;
- interf_props.Tref = TG;
- create_interface(dev, fluid, dummy_solid, &interf_props, &interf_TG);
-
- /* Create the TA interface */
- interf_props.temperature = SDIS_TEMPERATURE_NONE;
- interf_props.h = HA;
- interf_props.emissivity = 1;
- interf_props.Tref = TREF;
- create_interface(dev, solid, fluid_A, &interf_props, &interf_TA);
-
- /* Release the media */
- OK(sdis_medium_ref_put(solid));
- OK(sdis_medium_ref_put(dummy_solid));
- OK(sdis_medium_ref_put(fluid));
- OK(sdis_medium_ref_put(fluid_A));
-
- /* Front */
- model3d_interfaces[0] = interf_adiabatic_1;
- model3d_interfaces[1] = interf_adiabatic_1;
- model3d_interfaces[2] = interf_adiabatic_2;
- model3d_interfaces[3] = interf_adiabatic_2;
- /* Left */
- model3d_interfaces[4] = interf_TG;
- model3d_interfaces[5] = interf_TG;
- /* Back */
- model3d_interfaces[6] = interf_adiabatic_1;
- model3d_interfaces[7] = interf_adiabatic_1;
- model3d_interfaces[8] = interf_adiabatic_2;
- model3d_interfaces[9] = interf_adiabatic_2;
- /* Right */
- model3d_interfaces[10] = interf_TA;
- model3d_interfaces[11] = interf_TA;
- /* Top */
- model3d_interfaces[12] = interf_adiabatic_1;
- model3d_interfaces[13] = interf_adiabatic_1;
- model3d_interfaces[14] = interf_adiabatic_2;
- model3d_interfaces[15] = interf_adiabatic_2;
- /* Bottom */
- model3d_interfaces[16] = interf_adiabatic_1;
- model3d_interfaces[17] = interf_adiabatic_1;
- model3d_interfaces[18] = interf_adiabatic_2;
- model3d_interfaces[19] = interf_adiabatic_2;
- /* Inside */
- model3d_interfaces[20] = interf_P;
- model3d_interfaces[21] = interf_P;
-
- /* Bottom */
- model2d_interfaces[0] = interf_adiabatic_2;
- model2d_interfaces[1] = interf_adiabatic_1;
- /* Left */
- model2d_interfaces[2] = interf_TG;
- /* Top */
- model2d_interfaces[3] = interf_adiabatic_1;
- model2d_interfaces[4] = interf_adiabatic_2;
- /* Right */
- model2d_interfaces[5] = interf_TA;
- /* Contact */
- model2d_interfaces[6] = interf_P;
-
- /* Create the box scene */
- scn_args.get_indices = model3d_get_indices;
- scn_args.get_interface = model3d_get_interface;
- scn_args.get_position = model3d_get_position;
- scn_args.nprimitives = model3d_ntriangles;
- scn_args.nvertices = model3d_nvertices;
- scn_args.context = model3d_interfaces;
- scn_args.trad.temperature = TR;
- scn_args.trad.reference = TR;
- scn_args.t_range[0] = MMIN(MMIN(MMIN(MMIN(T0_FLUID, T0_SOLID), TA), TG), TR);
- scn_args.t_range[1] = MMAX(MMAX(MMAX(MMAX(T0_FLUID, T0_SOLID), TA), TG), TR);
- OK(sdis_scene_create(dev, &scn_args, &box_scn));
-
- /* Create the square scene */
- scn_args.get_indices = model2d_get_indices;
- scn_args.get_interface = model2d_get_interface;
- scn_args.get_position = model2d_get_position;
- scn_args.nprimitives = model2d_nsegments;
- scn_args.nvertices = model2d_nvertices;
- scn_args.context = model2d_interfaces;
- scn_args.trad.temperature = TR;
- scn_args.trad.reference = TR;
- scn_args.t_range[0] = MMIN(MMIN(MMIN(MMIN(T0_FLUID, T0_SOLID), TA), TG), TR);
- scn_args.t_range[1] = MMAX(MMAX(MMAX(MMAX(T0_FLUID, T0_SOLID), TA), TG), TR);
- OK(sdis_scene_2d_create(dev, &scn_args, &square_scn));
-
- /* Release the interfaces */
- OK(sdis_interface_ref_put(interf_adiabatic_1));
- OK(sdis_interface_ref_put(interf_adiabatic_2));
- OK(sdis_interface_ref_put(interf_TG));
- OK(sdis_interface_ref_put(interf_P));
- OK(sdis_interface_ref_put(interf_TA));
-
- /* Solve */
- OK(ssp_rng_create(NULL, SSP_RNG_KISS, &rng));
- printf(">> Box scene\n");
- solve_tfluid(box_scn);
- solve_tbound1(box_scn, rng);
- solve_tbound2(box_scn, rng);
- solve_tsolid(box_scn, rng);
- printf("\n>> Square scene\n");
- solve_tfluid(square_scn);
- solve_tbound1(box_scn, rng);
- solve_tbound2(box_scn, rng);
- solve_tsolid(square_scn, rng);
-
- OK(sdis_scene_ref_put(box_scn));
- OK(sdis_scene_ref_put(square_scn));
- OK(sdis_device_ref_put(dev));
- OK(ssp_rng_ref_put(rng));
-
- CHK(mem_allocated_size() == 0);
- return 0;
-}
-
diff --git a/src/test_sdis_unsteady_atm.c b/src/test_sdis_unsteady_atm.c
@@ -0,0 +1,956 @@
+/* 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 "sdis.h"
+#include "test_sdis_utils.h"
+
+#include <rsys/clock_time.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/double3.h>
+#include <rsys/math.h>
+#include <star/ssp.h>
+
+/*
+ * The physical configuration is the following: a slab of fluid with known
+ * thermophysical properties but unknown temperature is located between a
+ * "ground" and a slab of solid, with also a unknown temperature profile. On
+ * the other side of the solid slab, is an "atmosphere" with known temperature,
+ * and known radiative temperature.
+ *
+ * Solving the system means: finding the temperature of the ground, of the
+ * fluid, of the boundaries, and also the temperature inside the solid, at
+ * various locations (the 1D slab is discretized in order to obtain the
+ * reference)
+ *
+ * The reference for this system comes from a numerical method and is not
+ * analytic. Thus the compliance test MC VS reference is not the usual |MC -
+ * ref| <= 3*sigma but is |MC -ref| <= (Tmax -Tmin) * 0.01.
+ *
+ * 3D 2D
+ *
+ * /////////////// ///////////////
+ * +-----------+-+ +-----------+-+
+ * /' / /| | | |
+ * +-----------+-+ | HA _\ <---- TR | _\ | | HA _\ <---- TR
+ * | | _\ | | | / / <---- TR TG| HG / / HC | | / / <---- TR
+ * | |HG / / HC| | | TA \__/ <---- TR | \__/ | | TA \__/ <---- TR
+ * TG| | \__/ | | | | | |
+ * | +.........|.|.+ +-----------+-+
+ * |, |/|/ /////H///////E/
+ * +-----------+-+
+ * //////H//////E/
+ */
+
+#define XH 3
+#define XHpE 3.2
+#define XE (XHpE - XH)
+
+#define T0_SOLID 300
+#define T0_FLUID 300
+
+#define N 10000 /* #realisations */
+
+#define TG 310
+#define HG 400
+
+#define HC 400
+
+#define TA 290
+#define HA 400
+#define TR 260
+
+#define TMAX (MMAX(MMAX(MMAX(T0_FLUID, T0_SOLID), MMAX(TG, TA)), TR))
+#define TMIN (MMIN(MMIN(MMIN(T0_FLUID, T0_SOLID), MMIN(TG, TA)), TR))
+#define EPS ((TMAX-TMIN)*0.01)
+
+/* hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon
+ * Tref = (hr / (4 * 5.6696e-8 * epsilon)) ^ 1/3, hr = 6 */
+#define TREF 297.974852286
+
+#define RHO_F 1.3
+#define CP_F 1005
+#define RHO_S 2400
+#define CP_S 800
+#define LAMBDA 0.6
+
+#define X_PROBE (XH + 0.2 * XE)
+
+#define DELTA (XE/40.0)
+
+/*******************************************************************************
+ * Box geometry
+ ******************************************************************************/
+static const double model3d_vertices[12/*#vertices*/*3/*#coords per vertex*/] = {
+ 0, 0, 0,
+ XH, 0, 0,
+ XHpE, 0, 0,
+ 0, XHpE, 0,
+ XH, XHpE, 0,
+ XHpE, XHpE, 0,
+ 0, 0, XHpE,
+ XH, 0, XHpE,
+ XHpE, 0, XHpE,
+ 0, XHpE, XHpE,
+ XH, XHpE, XHpE,
+ XHpE, XHpE, XHpE
+};
+static const size_t model3d_nvertices = sizeof(model3d_vertices)/(sizeof(double)*3);
+
+/* The following array lists the indices toward the 3D vertices of each
+ * triangle.
+ * ,3---,4---,5 ,3----4----5 ,4
+ * ,' | ,' | ,'/| ,'/| \ | \ | ,'/|
+ * 9----10---11 / | 9' / | \ | \ | 10 / | Y
+ * |', |', | / ,2 | / ,0---,1---,2 | / ,1 |
+ * | ',| ',|/,' |/,' | ,' | ,' |/,' o--X
+ * 6----7----8' 6----7'---8' 7 /
+ * Front, right Back, left and Internal Z
+ * and Top faces bottom faces face */
+static const size_t model3d_indices[22/*#triangles*/*3/*#indices per triangle*/] = {
+ 0, 3, 1, 1, 3, 4, 1, 4, 2, 2, 4, 5, /* -Z */
+ 0, 6, 3, 3, 6, 9, /* -X */
+ 6, 7, 9, 9, 7, 10, 7, 8, 10, 10, 8, 11, /* +Z */
+ 5, 11, 8, 8, 2, 5, /* +X */
+ 3, 9, 10, 10, 4, 3, 4, 10, 11, 11, 5, 4, /* +Y */
+ 0, 1, 7, 7, 6, 0, 1, 2, 8, 8, 7, 1, /* -Y */
+ 4, 10, 7, 7, 1, 4 /* Inside */
+};
+static const size_t model3d_ntriangles = sizeof(model3d_indices)/(sizeof(size_t)*3);
+
+static INLINE void
+model3d_get_indices(const size_t itri, size_t ids[3], void* context)
+{
+ (void)context;
+ CHK(ids);
+ CHK(itri < model3d_ntriangles);
+ ids[0] = model3d_indices[itri * 3 + 0];
+ ids[1] = model3d_indices[itri * 3 + 1];
+ ids[2] = model3d_indices[itri * 3 + 2];
+}
+
+static INLINE void
+model3d_get_position(const size_t ivert, double pos[3], void* context)
+{
+ (void)context;
+ CHK(pos);
+ CHK(ivert < model3d_nvertices);
+ pos[0] = model3d_vertices[ivert * 3 + 0];
+ pos[1] = model3d_vertices[ivert * 3 + 1];
+ pos[2] = model3d_vertices[ivert * 3 + 2];
+}
+
+static INLINE void
+model3d_get_interface(const size_t itri, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ CHK(itri < model3d_ntriangles);
+ *bound = interfaces[itri];
+}
+
+/*******************************************************************************
+ * Square geometry
+ ******************************************************************************/
+static const double model2d_vertices[6/*#vertices*/ * 2/*#coords per vertex*/] = {
+ XHpE, 0,
+ XH, 0,
+ 0, 0,
+ 0, XHpE,
+ XH, XHpE,
+ XHpE, XHpE
+};
+static const size_t model2d_nvertices = sizeof(model2d_vertices)/(sizeof(double)*2);
+
+static const size_t model2d_indices[7/*#segments*/ * 2/*#indices per segment*/] = {
+ 0, 1, 1, 2, /* Bottom */
+ 2, 3, /* Left */
+ 3, 4, 4, 5, /* Top */
+ 5, 0, /* Right */
+ 4, 1 /* Inside */
+};
+static const size_t model2d_nsegments = sizeof(model2d_indices)/(sizeof(size_t)*2);
+
+static INLINE void
+model2d_get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ (void)context;
+ CHK(ids);
+ CHK(iseg < model2d_nsegments);
+ ids[0] = model2d_indices[iseg * 2 + 0];
+ ids[1] = model2d_indices[iseg * 2 + 1];
+}
+
+static INLINE void
+model2d_get_position(const size_t ivert, double pos[2], void* context)
+{
+ (void)context;
+ CHK(pos);
+ CHK(ivert < model2d_nvertices);
+ pos[0] = model2d_vertices[ivert * 2 + 0];
+ pos[1] = model2d_vertices[ivert * 2 + 1];
+}
+
+static INLINE void
+model2d_get_interface
+(const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ CHK(iseg < model2d_nsegments);
+ *bound = interfaces[iseg];
+}
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+struct solid {
+ double lambda;
+ double rho;
+ double cp;
+ double delta;
+ double temperature;
+ double t0;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ struct solid* solid;
+ CHK(vtx && data);
+ solid = ((struct solid*)sdis_data_cget(data));
+ return solid->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ struct solid* solid;
+ CHK(vtx && data);
+ solid = ((struct solid*)sdis_data_cget(data));
+ return solid->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ struct solid* solid;
+ CHK(vtx && data);
+ solid = ((struct solid*)sdis_data_cget(data));
+ return solid->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ struct solid* solid;
+ CHK(vtx && data);
+ solid = ((struct solid*)sdis_data_cget(data));
+ return solid->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ struct solid* solid;
+ CHK(vtx && data);
+ solid = ((struct solid*)sdis_data_cget(data));
+ if(vtx->time <= solid->t0)
+ return solid->temperature;
+ return SDIS_TEMPERATURE_NONE;
+}
+
+struct fluid {
+ double rho;
+ double cp;
+ double t0;
+ double temperature;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ struct fluid* fluid;
+ CHK(vtx && data);
+ fluid = ((struct fluid*)sdis_data_cget(data));
+ if(vtx->time <= fluid->t0)
+ return fluid->temperature;
+ return SDIS_TEMPERATURE_NONE;
+}
+
+static double
+fluid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ struct fluid* fluid;
+ CHK(vtx && data);
+ fluid = ((struct fluid*)sdis_data_cget(data));
+ return fluid->rho;
+}
+
+static double
+fluid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ struct fluid* fluid;
+ CHK(vtx && data);
+ fluid = ((struct fluid*)sdis_data_cget(data));
+ return fluid->cp;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double temperature;
+ double emissivity;
+ double h;
+ double Tref;
+};
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf;
+ CHK(frag && data);
+ interf = sdis_data_cget(data);
+ return interf->temperature;
+}
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf;
+ CHK(frag && data);
+ interf = sdis_data_cget(data);
+ return interf->h;
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
+{
+ const struct interf* interf;
+ (void)source_id;
+ CHK(frag && data);
+ interf = sdis_data_cget(data);
+ return interf->emissivity;
+}
+
+static double
+interface_get_Tref
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ const struct interf* interf;
+ CHK(frag && data);
+ interf = sdis_data_cget(data);
+ return interf->Tref;
+}
+
+/*******************************************************************************
+ * Radiative environment
+ ******************************************************************************/
+static double
+radenv_get_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray, (void)data;
+ return TR; /* [K] */
+}
+
+static double
+radenv_get_reference_temperature
+ (const struct sdis_radiative_ray* ray,
+ struct sdis_data* data)
+{
+ (void)ray, (void)data;
+ return TR; /* [K] */
+}
+
+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;
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+create_interface
+ (struct sdis_device* dev,
+ struct sdis_medium* front,
+ struct sdis_medium* back,
+ const struct interf* interf,
+ struct sdis_interface** out_interf)
+{
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_data* data = NULL;
+
+ CHK(interf != NULL);
+
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) {
+ shader.convection_coef = interface_get_convection_coef;
+ shader.convection_coef_upper_bound = interf->h;
+ }
+ if(sdis_medium_get_type(front) == SDIS_FLUID) {
+ shader.front.emissivity = interface_get_emissivity;
+ shader.front.reference_temperature = interface_get_Tref;
+ }
+ if(sdis_medium_get_type(back) == SDIS_FLUID) {
+ shader.back.emissivity = interface_get_emissivity;
+ shader.back.reference_temperature = interface_get_Tref;
+ }
+
+ OK(sdis_data_create(dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data));
+ *((struct interf*)sdis_data_get(data)) = *interf;
+
+ OK(sdis_interface_create(dev, front, back, &shader, data, out_interf));
+ OK(sdis_data_ref_put(data));
+}
+
+static void
+solve_tbound1
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng)
+{
+ char dump[128];
+ struct time t0, t1;
+ struct sdis_estimator* estimator;
+ struct sdis_solve_probe_boundary_args solve_args
+ = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+ enum sdis_scene_dimension dim;
+ const double t[] = { 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
+ const double ref[sizeof(t) / sizeof(*t)] = {
+ 290.046375, 289.903935, 289.840490, 289.802690, 289.777215, 289.759034,
+ 289.745710, 289.735826, 289.728448, 289.722921
+ };
+ const int nsimuls = sizeof(t) / sizeof(*t);
+ int isimul;
+ ASSERT(scn && rng);
+
+ OK(sdis_scene_get_dimension(scn, &dim));
+
+ solve_args.nrealisations = N;
+ solve_args.side = SDIS_FRONT;
+ FOR_EACH(isimul, 0, nsimuls) {
+ solve_args.time_range[0] = solve_args.time_range[1] = t[isimul];
+ if(dim == SDIS_SCENE_2D) {
+ solve_args.iprim = model2d_nsegments - 1;
+ solve_args.uv[0] = ssp_rng_uniform_double(rng, 0, 1);
+ } else {
+ double u = ssp_rng_uniform_double(rng, 0, 1);
+ double v = ssp_rng_uniform_double(rng, 0, 1);
+ double w = ssp_rng_uniform_double(rng, 0, 1);
+ double x = 1 / (u + v + w);
+ solve_args.iprim = (isimul % 2) ? 10 : 11; /* +X face */
+ solve_args.uv[0] = u * x;
+ solve_args.uv[1] = v * x;
+ }
+
+ time_current(&t0);
+ OK(sdis_solve_probe_boundary(scn, &solve_args, &estimator));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
+ switch(dim) {
+ case SDIS_SCENE_2D:
+ printf("Unstationary temperature at (%lu/%g) at t=%g = %g ~ %g +/- %g\n",
+ (unsigned long)solve_args.iprim, solve_args.uv[0], t[isimul],
+ ref[isimul], T.E, T.SE);
+ break;
+ case SDIS_SCENE_3D:
+ printf("Unstationary temperature at (%lu/%g,%g) at t=%g = %g ~ %g +/- %g\n",
+ (unsigned long)solve_args.iprim, SPLIT2(solve_args.uv), t[isimul],
+ ref[isimul], T.E, T.SE);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ printf("Elapsed time = %s\n", dump);
+ printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
+
+ CHK(eq_eps(T.E, ref[isimul], EPS));
+ /*CHK(eq_eps(T.E, ref[isimul], T.SE*3));*/
+
+ OK(sdis_estimator_ref_put(estimator));
+ }
+}
+
+static void
+solve_tbound2
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng)
+{
+ char dump[128];
+ struct time t0, t1;
+ struct sdis_estimator* estimator;
+ struct sdis_solve_probe_boundary_args solve_args
+ = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+ enum sdis_scene_dimension dim;
+ const double t[] = { 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
+ const double ref[sizeof(t) / sizeof(*t)] = {
+ 309.08032, 309.34626, 309.46525, 309.53625, 309.58408, 309.618121,
+ 309.642928, 309.661167, 309.674614, 309.684524
+ };
+ const int nsimuls = sizeof(t) / sizeof(*t);
+ int isimul;
+ ASSERT(scn && rng);
+
+ OK(sdis_scene_get_dimension(scn, &dim));
+
+ solve_args.nrealisations = N;
+ solve_args.side = SDIS_FRONT;
+ FOR_EACH(isimul, 0, nsimuls) {
+ solve_args.time_range[0] = solve_args.time_range[1] = t[isimul];
+ if(dim == SDIS_SCENE_2D) {
+ solve_args.iprim = model2d_nsegments - 1;
+ solve_args.uv[0] = ssp_rng_uniform_double(rng, 0, 1);
+ } else {
+ double u = ssp_rng_uniform_double(rng, 0, 1);
+ double v = ssp_rng_uniform_double(rng, 0, 1);
+ double w = ssp_rng_uniform_double(rng, 0, 1);
+ double x = 1 / (u + v + w);
+ solve_args.iprim = (isimul % 2) ? 20 : 21; /* Internal face */
+ solve_args.uv[0] = u * x;
+ solve_args.uv[1] = v * x;
+ }
+
+ time_current(&t0);
+ OK(sdis_solve_probe_boundary(scn, &solve_args, &estimator));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
+ switch(dim) {
+ case SDIS_SCENE_2D:
+ printf("Unstationary temperature at (%lu/%g) at t=%g = %g ~ %g +/- %g\n",
+ (unsigned long)solve_args.iprim, solve_args.uv[0], t[isimul],
+ ref[isimul], T.E, T.SE);
+ break;
+ case SDIS_SCENE_3D:
+ printf("Unstationary temperature at (%lu/%g,%g) at t=%g = %g ~ %g +/- %g\n",
+ (unsigned long)solve_args.iprim, SPLIT2(solve_args.uv), t[isimul],
+ ref[isimul], T.E, T.SE);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ printf("Elapsed time = %s\n", dump);
+ printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails <= N/1000);
+ CHK(eq_eps(T.E, ref[isimul], EPS));
+ /*CHK(eq_eps(T.E, ref[isimul], T.SE*3));*/
+
+ OK(sdis_estimator_ref_put(estimator));
+ }
+}
+
+static void
+solve_tsolid
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng)
+{
+ char dump[128];
+ struct time t0, t1;
+ struct sdis_estimator* estimator;
+ struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+ enum sdis_scene_dimension dim;
+ const double t[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
+ const double ref[sizeof(t) / sizeof(*t)] = {
+ 300, 300.87408, 302.25832, 303.22164, 303.89954, 304.39030, 304.75041,
+ 305.01595, 305.21193, 305.35641, 305.46271
+ };
+ const int nsimuls = sizeof(t) / sizeof(*t);
+ int isimul;
+ ASSERT(scn && rng);
+
+ OK(sdis_scene_get_dimension(scn, &dim));
+
+ solve_args.nrealisations = N;
+ FOR_EACH(isimul, 0, nsimuls) {
+ solve_args.time_range[0] = solve_args.time_range[1] = t[isimul];
+ solve_args.position[0] = X_PROBE;
+ solve_args.position[1] = ssp_rng_uniform_double(rng, 0.1*XHpE, 0.9*XHpE);
+
+ if(dim == SDIS_SCENE_3D)
+ solve_args.position[2] = ssp_rng_uniform_double(rng, 0.1*XHpE, 0.9*XHpE);
+
+ time_current(&t0);
+ OK(sdis_solve_probe(scn, &solve_args, &estimator));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
+ switch (dim) {
+ case SDIS_SCENE_2D:
+ printf("Unstationary temperature at (%g,%g) at t=%g = %g ~ %g +/- %g\n",
+ SPLIT2(solve_args.position), t[isimul], ref[isimul], T.E, T.SE);
+ break;
+ case SDIS_SCENE_3D:
+ printf("Unstationary temperature at (%g,%g,%g) at t=%g = %g ~ %g +/- %g\n",
+ SPLIT3(solve_args.position), t[isimul], ref[isimul], T.E, T.SE);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ printf("Elapsed time = %s\n", dump);
+ printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails <= N / 1000);
+ CHK(eq_eps(T.E, ref[isimul], EPS));
+ /*CHK(eq_eps(T.E, ref[isimul], T.SE*3));*/
+
+ OK(sdis_estimator_ref_put(estimator));
+ }
+}
+
+static void
+solve_tfluid
+ (struct sdis_scene* scn)
+{
+ char dump[128];
+ struct time t0, t1;
+ struct sdis_estimator* estimator;
+ struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+ enum sdis_scene_dimension dim;
+ double eps;
+ const double t[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 };
+ const double ref[sizeof(t) / sizeof(*t)] = {
+ 300, 309.53905, 309.67273, 309.73241, 309.76798, 309.79194, 309.80899,
+ 309.82141, 309.83055, 309.83728, 309.84224
+ };
+ const int nsimuls = sizeof(t) / sizeof(*t);
+ int isimul;
+ ASSERT(scn);
+
+ OK(sdis_scene_get_dimension(scn, &dim));
+
+ solve_args.nrealisations = N;
+ solve_args.position[0] = XH * 0.5;
+ solve_args.position[1] = XH * 0.5;
+ solve_args.position[2] = XH * 0.5;
+ FOR_EACH(isimul, 0, nsimuls) {
+ solve_args.time_range[0] = solve_args.time_range[1] = t[isimul];
+
+ time_current(&t0);
+ OK(sdis_solve_probe(scn, &solve_args, &estimator));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
+ switch (dim) {
+ case SDIS_SCENE_2D:
+ printf("Unstationary fluid temperature at t=%g = %g ~ %g +/- %g\n",
+ t[isimul], ref[isimul], T.E, T.SE);
+ break;
+ case SDIS_SCENE_3D:
+ printf("Unstationary fluid temperature at t=%g = %g ~ %g +/- %g\n",
+ t[isimul], ref[isimul], T.E, T.SE);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ printf("Elapsed time = %s\n", dump);
+ printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails <= N / 1000);
+
+ eps = EPS;
+ CHK(eq_eps(T.E, ref[isimul], eps));
+
+ OK(sdis_estimator_ref_put(estimator));
+ }
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct sdis_data* data = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* fluid_A = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* dummy_solid = NULL;
+ struct sdis_interface* interf_adiabatic_1 = NULL;
+ struct sdis_interface* interf_adiabatic_2 = NULL;
+ struct sdis_interface* interf_TG = NULL;
+ struct sdis_interface* interf_P = NULL;
+ struct sdis_interface* interf_TA = NULL;
+ struct sdis_radiative_env* radenv = NULL;
+ struct sdis_scene* box_scn = NULL;
+ struct sdis_scene* square_scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface* model3d_interfaces[22 /*#triangles*/];
+ struct sdis_interface* model2d_interfaces[7/*#segments*/];
+ struct interf interf_props;
+ struct solid* solid_props = NULL;
+ struct fluid* fluid_props = NULL;
+ struct ssp_rng* rng = NULL;
+ (void)argc, (void)argv;
+
+ OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev));
+
+ radenv = create_radenv(dev);
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+
+ /* Create the solid media */
+ OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data));
+ solid_props = sdis_data_get(data);
+ solid_props->lambda = LAMBDA;
+ solid_props->cp = CP_S;
+ solid_props->rho = RHO_S;
+ solid_props->delta = DELTA;
+ solid_props->t0 = 0;
+ solid_props->temperature = T0_SOLID;
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid));
+ OK(sdis_data_ref_put(data));
+
+ /* Create a dummy solid media to be used outside the model */
+ OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data));
+ solid_props = sdis_data_get(data);
+ solid_props->lambda = 0;
+ solid_props->cp = 1;
+ solid_props->rho = 1;
+ solid_props->delta = 1;
+ solid_props->t0 = INF;
+ solid_props->temperature = SDIS_TEMPERATURE_NONE;
+ OK(sdis_solid_create(dev, &solid_shader, data, &dummy_solid));
+ OK(sdis_data_ref_put(data));
+
+ /* Setup the fluid shader */
+ fluid_shader.calorific_capacity = fluid_get_calorific_capacity;
+ fluid_shader.volumic_mass = fluid_get_volumic_mass;
+ fluid_shader.temperature = fluid_get_temperature;
+
+ /* Create the internal fluid media */
+ OK(sdis_data_create(dev, sizeof(struct fluid), 16, NULL, &data));
+ fluid_props = sdis_data_get(data);
+ fluid_props->cp = CP_F;
+ fluid_props->rho = RHO_F;
+ fluid_props->t0 = 0;
+ fluid_props->temperature = T0_FLUID;
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the 'A' fluid media */
+ OK(sdis_data_create(dev, sizeof(struct fluid), 16, NULL, &data));
+ fluid_props = sdis_data_get(data);
+ fluid_props->cp = 1;
+ fluid_props->rho = 1;
+ fluid_props->t0 = INF;
+ fluid_props->temperature = TA;
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid_A));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the adiabatic interfaces */
+ interf_props.temperature = SDIS_TEMPERATURE_NONE;
+ interf_props.h = 0;
+ interf_props.emissivity = 0;
+ interf_props.Tref = TREF;
+ create_interface(dev, fluid, dummy_solid, &interf_props, &interf_adiabatic_1);
+ create_interface(dev, solid, dummy_solid, &interf_props, &interf_adiabatic_2);
+
+ /* Create the P interface */
+ interf_props.temperature = SDIS_TEMPERATURE_NONE;
+ interf_props.h = HC;
+ interf_props.emissivity = 1;
+ interf_props.Tref = TREF;
+ create_interface(dev, fluid, solid, &interf_props, &interf_P);
+
+ /* Create the TG interface */
+ interf_props.temperature = TG;
+ interf_props.h = HG;
+ interf_props.emissivity = 1;
+ interf_props.Tref = TG;
+ create_interface(dev, fluid, dummy_solid, &interf_props, &interf_TG);
+
+ /* Create the TA interface */
+ interf_props.temperature = SDIS_TEMPERATURE_NONE;
+ interf_props.h = HA;
+ interf_props.emissivity = 1;
+ interf_props.Tref = TREF;
+ create_interface(dev, solid, fluid_A, &interf_props, &interf_TA);
+
+ /* Release the media */
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_medium_ref_put(dummy_solid));
+ OK(sdis_medium_ref_put(fluid));
+ OK(sdis_medium_ref_put(fluid_A));
+
+ /* Front */
+ model3d_interfaces[0] = interf_adiabatic_1;
+ model3d_interfaces[1] = interf_adiabatic_1;
+ model3d_interfaces[2] = interf_adiabatic_2;
+ model3d_interfaces[3] = interf_adiabatic_2;
+ /* Left */
+ model3d_interfaces[4] = interf_TG;
+ model3d_interfaces[5] = interf_TG;
+ /* Back */
+ model3d_interfaces[6] = interf_adiabatic_1;
+ model3d_interfaces[7] = interf_adiabatic_1;
+ model3d_interfaces[8] = interf_adiabatic_2;
+ model3d_interfaces[9] = interf_adiabatic_2;
+ /* Right */
+ model3d_interfaces[10] = interf_TA;
+ model3d_interfaces[11] = interf_TA;
+ /* Top */
+ model3d_interfaces[12] = interf_adiabatic_1;
+ model3d_interfaces[13] = interf_adiabatic_1;
+ model3d_interfaces[14] = interf_adiabatic_2;
+ model3d_interfaces[15] = interf_adiabatic_2;
+ /* Bottom */
+ model3d_interfaces[16] = interf_adiabatic_1;
+ model3d_interfaces[17] = interf_adiabatic_1;
+ model3d_interfaces[18] = interf_adiabatic_2;
+ model3d_interfaces[19] = interf_adiabatic_2;
+ /* Inside */
+ model3d_interfaces[20] = interf_P;
+ model3d_interfaces[21] = interf_P;
+
+ /* Bottom */
+ model2d_interfaces[0] = interf_adiabatic_2;
+ model2d_interfaces[1] = interf_adiabatic_1;
+ /* Left */
+ model2d_interfaces[2] = interf_TG;
+ /* Top */
+ model2d_interfaces[3] = interf_adiabatic_1;
+ model2d_interfaces[4] = interf_adiabatic_2;
+ /* Right */
+ model2d_interfaces[5] = interf_TA;
+ /* Contact */
+ model2d_interfaces[6] = interf_P;
+
+ /* Create the box scene */
+ scn_args.get_indices = model3d_get_indices;
+ scn_args.get_interface = model3d_get_interface;
+ scn_args.get_position = model3d_get_position;
+ scn_args.nprimitives = model3d_ntriangles;
+ scn_args.nvertices = model3d_nvertices;
+ scn_args.context = model3d_interfaces;
+ scn_args.radenv = radenv;
+ scn_args.t_range[0] = MMIN(MMIN(MMIN(MMIN(T0_FLUID, T0_SOLID), TA), TG), TR);
+ scn_args.t_range[1] = MMAX(MMAX(MMAX(MMAX(T0_FLUID, T0_SOLID), TA), TG), TR);
+ OK(sdis_scene_create(dev, &scn_args, &box_scn));
+
+ /* Create the square scene */
+ scn_args.get_indices = model2d_get_indices;
+ scn_args.get_interface = model2d_get_interface;
+ scn_args.get_position = model2d_get_position;
+ scn_args.nprimitives = model2d_nsegments;
+ scn_args.nvertices = model2d_nvertices;
+ scn_args.context = model2d_interfaces;
+ scn_args.radenv = radenv;
+ scn_args.t_range[0] = MMIN(MMIN(MMIN(MMIN(T0_FLUID, T0_SOLID), TA), TG), TR);
+ scn_args.t_range[1] = MMAX(MMAX(MMAX(MMAX(T0_FLUID, T0_SOLID), TA), TG), TR);
+ OK(sdis_scene_2d_create(dev, &scn_args, &square_scn));
+
+ /* Release the interfaces */
+ OK(sdis_interface_ref_put(interf_adiabatic_1));
+ OK(sdis_interface_ref_put(interf_adiabatic_2));
+ OK(sdis_interface_ref_put(interf_TG));
+ OK(sdis_interface_ref_put(interf_P));
+ OK(sdis_interface_ref_put(interf_TA));
+
+ /* Solve */
+ OK(ssp_rng_create(NULL, SSP_RNG_KISS, &rng));
+ printf(">> Box scene\n");
+ solve_tfluid(box_scn);
+ solve_tbound1(box_scn, rng);
+ solve_tbound2(box_scn, rng);
+ solve_tsolid(box_scn, rng);
+ printf("\n>> Square scene\n");
+ solve_tfluid(square_scn);
+ solve_tbound1(box_scn, rng);
+ solve_tbound2(box_scn, rng);
+ solve_tsolid(square_scn, rng);
+
+ OK(sdis_radiative_env_ref_put(radenv));
+ OK(sdis_scene_ref_put(box_scn));
+ OK(sdis_scene_ref_put(square_scn));
+ OK(sdis_device_ref_put(dev));
+ OK(ssp_rng_ref_put(rng));
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -86,19 +86,23 @@ accum_flux_terms
static res_T
solve_green_path(struct sdis_green_path* path, void* ctx)
{
- struct sdis_point pt = SDIS_POINT_NULL;
+ struct sdis_green_path_end end = SDIS_GREEN_PATH_END_NULL;
+
struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+ struct sdis_radiative_ray ray = SDIS_RADIATIVE_RAY_NULL;
+
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_radiative_env_shader radenv = SDIS_RADIATIVE_ENV_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 power = 0;
double flux = 0;
double external_flux = 0; /* [W/m^2] */
@@ -143,48 +147,41 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
external_flux *= sdis_source_get_power(source, INF); /* [W] */
}
- 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));
- OK(sdis_green_path_get_end_type(path, &end_type));
-
- BA(sdis_green_path_get_limit_point(NULL, NULL));
- BA(sdis_green_path_get_limit_point(NULL, &pt));
- BA(sdis_green_path_get_limit_point(path, NULL));
- if(end_type == SDIS_GREEN_PATH_END_RADIATIVE) {
- struct sdis_ambient_radiative_temperature trad;
- BO(sdis_green_path_get_limit_point(path, &pt));
-
- BA(sdis_scene_get_ambient_radiative_temperature(NULL, NULL));
- 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));
- temp = trad.temperature;
- } else {
- OK(sdis_green_path_get_limit_point(path, &pt));
- switch(pt.type) {
- case SDIS_FRAGMENT:
- frag = pt.data.itfrag.fragment;
- OK(sdis_interface_get_shader(pt.data.itfrag.intface, &interf));
- data = sdis_interface_get_data(pt.data.itfrag.intface);
+ BA(sdis_green_path_get_end(NULL, NULL));
+ BA(sdis_green_path_get_end(NULL, &end));
+ BA(sdis_green_path_get_end(path, NULL));
+ OK(sdis_green_path_get_end(path, &end));
+ switch(end.type) {
+ case SDIS_GREEN_PATH_END_AT_INTERFACE:
+ frag = end.data.itfrag.fragment;
+ OK(sdis_interface_get_shader(end.data.itfrag.intface, &interf));
+ data = sdis_interface_get_data(end.data.itfrag.intface);
temp = frag.side == SDIS_FRONT
? interf.front.temperature(&frag, data)
: interf.back.temperature(&frag, data);
break;
- case SDIS_VERTEX:
- vtx = pt.data.mdmvert.vertex;
- type = sdis_medium_get_type(pt.data.mdmvert.medium);
- data = sdis_medium_get_data(pt.data.mdmvert.medium);
+
+ case SDIS_GREEN_PATH_END_AT_RADIATIVE_ENV:
+ ray = end.data.radenvray.ray;
+ OK(sdis_radiative_env_get_shader(end.data.radenvray.radenv, &radenv));
+ data = sdis_radiative_env_get_data(end.data.radenvray.radenv);
+ temp = radenv.temperature(&ray, data);
+ break;
+
+ case SDIS_GREEN_PATH_END_IN_VOLUME:
+ vtx = end.data.mdmvert.vertex;
+ type = sdis_medium_get_type(end.data.mdmvert.medium);
+ data = sdis_medium_get_data(end.data.mdmvert.medium);
if(type == SDIS_FLUID) {
- OK(sdis_fluid_get_shader(pt.data.mdmvert.medium, &fluid));
+ OK(sdis_fluid_get_shader(end.data.mdmvert.medium, &fluid));
temp = fluid.temperature(&vtx, data);
} else {
- OK(sdis_solid_get_shader(pt.data.mdmvert.medium, &solid));
+ OK(sdis_solid_get_shader(end.data.mdmvert.medium, &solid));
temp = solid.temperature(&vtx, data);
}
break;
+
default: FATAL("Unreachable code.\n"); break;
- }
}
weight = temp + power + external_flux + flux;
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -149,7 +149,7 @@ square_get_interface
}
/*******************************************************************************
- * Medium & interface
+ * Medium, interface and ray
******************************************************************************/
static INLINE double
dummy_medium_getter
@@ -169,6 +169,25 @@ dummy_interface_getter
return 0;
}
+static INLINE double
+dummy_radiative_interface_getter
+ (const struct sdis_interface_fragment* frag,
+ const unsigned source_id,
+ struct sdis_data* data)
+{
+ (void)data, (void)source_id;
+ CHK(frag != NULL);
+ return 0;
+}
+
+static INLINE double
+dummy_ray_getter(const struct sdis_radiative_ray* ray, struct sdis_data* data)
+{
+ (void)data;
+ CHK(ray != NULL);
+ return 0;
+}
+
static const struct sdis_solid_shader DUMMY_SOLID_SHADER = {
dummy_medium_getter, /* Calorific capacity */
dummy_medium_getter, /* Thermal conductivity */
@@ -186,12 +205,11 @@ static const struct sdis_fluid_shader DUMMY_FLUID_SHADER = {
0 /* Initial time */
};
-
#define DUMMY_INTERFACE_SIDE_SHADER__ { \
dummy_interface_getter, /* Temperature */ \
dummy_interface_getter, /* Flux */ \
- dummy_interface_getter, /* Emissivity */ \
- dummy_interface_getter, /* Specular fraction */ \
+ dummy_radiative_interface_getter, /* Emissivity */ \
+ dummy_radiative_interface_getter, /* Specular fraction */ \
dummy_interface_getter, /* Reference temperature */ \
1 /* Handle external flux */ \
}
@@ -203,6 +221,11 @@ static const struct sdis_interface_shader DUMMY_INTERFACE_SHADER = {
DUMMY_INTERFACE_SIDE_SHADER__ /* Back side */
};
+static const struct sdis_radiative_env_shader DUMMY_RAY_SHADER = {
+ dummy_ray_getter,
+ dummy_ray_getter
+};
+
/*******************************************************************************
* Device creation
******************************************************************************/
