commit 436851de4649b77ccca27d8075222dda967c42e0
parent 252fdf5a9c1e14aadbd66df847f60e91398a13a9
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 10 Apr 2024 18:24:08 +0200
Improved definition of the radiative environment
This commit is a complete rewrite of the way the radiative environment
is managed. Until now, it was defined by two constants, the ambient
radiative temperature and a reference temperature used during
linearization of radiative transfer. These constants were then given as
input arguments during scene creation.
From now on, the radiative environment is described by its own
API-opaque data structure, as are interfaces, media and external
sources. A pointer to a variable in the radiative environment has
replaced the previous constants as the input argument to scene creation.
In addition, the temperature and reference of the radiative environment
are no longer constants. They are retrieved by functions dependent on a
given direction and source.
In fact, in addition to new functionalities, this profound change
simplifies the way in which the radiative environment is defined and
processed. It now follows the same API as other data (media, interfaces,
sources). For example, the green function doesn't have to handle a
specific case when a path ends in the radiative environment. It is now
handled like any other path boundary, as when it ends in a media or at
an interface.
Diffstat:
24 files changed, 1829 insertions(+), 1283 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\
@@ -197,11 +198,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 =\
@@ -303,10 +304,10 @@ 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 \
@@ -337,10 +338,10 @@ 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 \
@@ -371,10 +372,10 @@ 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,15 @@ 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 */
+ unsigned source_id; /* Identifier of the source */
+};
+#define SDIS_RADIATIVE_RAY_NULL__ {{0,0,0}, SDIS_INTERN_SOURCE_ID}
+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 +249,12 @@ typedef double
(const struct sdis_interface_fragment* frag, /* Interface position */
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 */
@@ -326,6 +347,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 +400,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 +483,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 +495,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 +502,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 +516,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 +1066,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.
******************************************************************************/
@@ -1141,19 +1193,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
@@ -1259,6 +1298,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
******************************************************************************/
@@ -1409,18 +1454,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_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_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;
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 */
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(radenv && ray && d3_is_normalized(ray->dir));
+ return 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(radenv && ray && d3_is_normalized(ray->dir));
+ return 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/test_sdis_draw_external_flux.c b/src/test_sdis_draw_external_flux.c
@@ -204,6 +204,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 +278,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 +293,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 +375,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 +391,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 +410,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
@@ -283,6 +283,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 +383,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 +399,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.source = source;
+ scn_args.radenv = radenv;
scn_args.context = &context;
OK(sdis_scene_2d_create(sdis, &scn_args, &scn));
@@ -386,7 +419,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 +435,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 +524,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 +542,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 +558,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
@@ -347,19 +347,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 +422,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 +463,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 +499,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 +553,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_picard.c b/src/test_sdis_picard.c
@@ -370,6 +370,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 +528,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 +563,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 +572,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 +603,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 +617,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 +639,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 +691,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 +704,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 +718,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 +732,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 +753,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 +779,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 +793,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_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
@@ -182,6 +182,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 +268,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 +284,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 +294,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 +381,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 +393,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 +609,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 +617,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
@@ -284,6 +284,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 +342,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 +351,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 +377,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 +386,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 +411,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 +422,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 +432,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 +604,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 +626,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 +686,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 +723,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 +813,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 +825,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_probe.c b/src/test_sdis_solve_probe.c
@@ -179,6 +179,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 +323,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 +339,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 +364,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 +373,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 +390,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 +408,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 +421,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 +496,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 +511,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 +538,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 +548,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 +625,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 +649,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 +669,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 +683,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_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,952 @@
+/* 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;
+}
+
+/*******************************************************************************
+ * 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,14 @@ dummy_interface_getter
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,7 +194,6 @@ 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 */ \
@@ -203,6 +210,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
******************************************************************************/
