commit bb0c8f4c24b48cd1630feffe15b8e5940c69930b
parent 74edf309a4297dcdec319c131908e7eba2039fc6
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 17 Sep 2020 14:01:56 +0200
Merge branch 'release_0.10'
Diffstat:
44 files changed, 2723 insertions(+), 166 deletions(-)
diff --git a/README.md b/README.md
@@ -25,7 +25,22 @@ variable the install directories of its dependencies.
## Release notes
-### Version 0.9.0
+### Version 0.10
+
+- Add support of green function [de]serialization. The
+ `sdis_green_function_write` function serializes the green function into a
+ stream while the `sdis_green_function_create_from_stream` function
+ deserialize it. Note that the scene used to deserialize the green function
+ must be the same of the one used to estimate it: the media and the interfaces
+ have to be created in the same order, the scene geometry must be the same,
+ etc.
+- Add the `sdis_scene_find_closest_point` function: search the point onto the
+ scene geometry that is the closest of the submitted position.
+- Add the `sdis_compute_power` function that evaluates the power of a medium.
+- Update the solver: the time of the sampled path is now rewind on solid
+ reinjection.
+
+### Version 0.9
- Update the API of the solve functions: the parameters of the simulation are
now grouped into a unique data structure rather than separately submitted as
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -29,12 +29,12 @@ CMAKE_DEPENDENT_OPTION(ALL_TESTS
# Check dependencies
###############################################################################
find_package(RCMake 0.4 REQUIRED)
-find_package(Star2D 0.3.1 REQUIRED)
-find_package(Star3D 0.6.2 REQUIRED)
+find_package(Star2D 0.4 REQUIRED)
+find_package(Star3D 0.7 REQUIRED)
find_package(StarSP 0.8 REQUIRED)
find_package(StarEnc2D 0.5 REQUIRED)
find_package(StarEnc3D 0.5 REQUIRED)
-find_package(RSys 0.8.1 REQUIRED)
+find_package(RSys 0.10 REQUIRED)
find_package(OpenMP 2.0 REQUIRED)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
@@ -56,7 +56,7 @@ rcmake_append_runtime_dirs(_runtime_dirs
# Configure and define targets
###############################################################################
set(VERSION_MAJOR 0)
-set(VERSION_MINOR 9)
+set(VERSION_MINOR 10)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
@@ -71,6 +71,7 @@ set(SDIS_FILES_SRC
sdis_interface.c
sdis_log.c
sdis_medium.c
+ sdis_misc.c
sdis_realisation.c
sdis_scene.c
sdis_solve.c)
@@ -90,6 +91,8 @@ set(SDIS_FILES_INC
sdis_heat_path_radiative_Xd.h
sdis_interface_c.h
sdis_log.h
+ sdis_misc.h
+ sdis_misc_Xd.h
sdis_medium_c.h
sdis_realisation.h
sdis_realisation_Xd.h
@@ -172,6 +175,7 @@ if(NOT NO_TEST)
endfunction()
new_test(test_sdis_camera)
+ new_test(test_sdis_compute_mean_power)
new_test(test_sdis_conducto_radiative)
new_test(test_sdis_conducto_radiative_2d)
new_test(test_sdis_convection)
@@ -194,6 +198,7 @@ if(NOT NO_TEST)
new_test(test_sdis_solve_boundary_flux)
new_test(test_sdis_solve_medium)
new_test(test_sdis_solve_medium_2d)
+ new_test(test_sdis_transcient)
new_test(test_sdis_volumic_power)
new_test(test_sdis_volumic_power4)
@@ -208,6 +213,7 @@ if(NOT NO_TEST)
add_test(test_sdis_volumic_power3_2d test_sdis_volumic_power3_2d)
endif()
+ target_link_libraries(test_sdis_compute_mean_power Star3DUT)
target_link_libraries(test_sdis_solid_random_walk_robustness Star3DUT)
target_link_libraries(test_sdis_solve_medium Star3DUT)
target_link_libraries(test_sdis_solve_probe3 Star3DUT)
diff --git a/src/sdis.h b/src/sdis.h
@@ -43,6 +43,7 @@
#define SDIS_VOLUMIC_POWER_NONE 0 /* <=> No volumic power */
#define SDIS_FLUX_NONE DBL_MAX /* <=> No flux */
+#define SDIS_PRIMITIVE_NONE SIZE_MAX /* Invalid primitive */
/* Forward declaration of external opaque data types */
struct logger;
@@ -113,8 +114,9 @@ static const struct sdis_interface_fragment SDIS_INTERFACE_FRAGMENT_NULL =
* Estimation data types
******************************************************************************/
enum sdis_estimator_type {
- SDIS_ESTIMATOR_TEMPERATURE,
- SDIS_ESTIMATOR_FLUX,
+ SDIS_ESTIMATOR_TEMPERATURE, /* In Kelvin */
+ SDIS_ESTIMATOR_FLUX, /* In Watt/m^2 */
+ SDIS_ESTIMATOR_POWER, /* In Watt */
SDIS_ESTIMATOR_TYPES_COUNT__
};
@@ -501,6 +503,23 @@ struct sdis_solve_camera_args {
static const struct sdis_solve_camera_args SDIS_SOLVE_CAMERA_ARGS_DEFAULT =
SDIS_SOLVE_CAMERA_ARGS_DEFAULT__;
+struct sdis_compute_power_args {
+ size_t nrealisations;
+ struct sdis_medium* medium; /* Medium to solve */
+ double time_range[2]; /* Observation time */
+ double fp_to_meter; /* Scale from floating point units to meters */
+ struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
+};
+#define SDIS_COMPUTE_POWER_ARGS_DEFAULT__ { \
+ 10000, /* #realisations */ \
+ NULL, /* Medium */ \
+ {DBL_MAX,DBL_MAX}, /* Time range */ \
+ 1, /* FP to meter */ \
+ NULL /* RNG state */ \
+}
+static const struct sdis_compute_power_args
+SDIS_COMPUTE_POWER_ARGS_DEFAULT = SDIS_COMPUTE_POWER_ARGS_DEFAULT__;
+
BEGIN_DECLS
/*******************************************************************************
@@ -791,6 +810,20 @@ sdis_scene_get_aabb
double lower[3],
double upper[3]);
+/* Search the point onto the scene geometry that is the closest of `pos'. The
+ * `radius' parameter controls the maximum search distance around `pos'. The
+ * returned closest point is expressed locally to the geometric primitive onto
+ * which it lies. If not found, the returned primitive is SDIS_PRIMITIVE_NONE.
+ * Note that even though only one point is returned, several position can have
+ * the same minimal distance to the queried position. */
+SDIS_API res_T
+sdis_scene_find_closest_point
+ (const struct sdis_scene* scn,
+ const double pos[3], /* Query position */
+ const double radius, /* Maximum search distance around pos */
+ size_t* iprim, /* Primitive index onto which the closest point lies */
+ double uv[2]); /* Parametric cordinate onto the primitive */
+
/* Define the world space position of a point onto the primitive `iprim' whose
* parametric coordinate is uv. */
SDIS_API res_T
@@ -800,7 +833,7 @@ sdis_scene_get_boundary_position
const double uv[2], /* Parametric coordinate onto the primitive */
double pos[3]); /* World space position */
-/* Project a world space position onto a primitive wrt its normal and compute
+/* roject a world space position onto a primitive wrt its normal and compute
* the parametric coordinates of the projected point onto the primitive. This
* function may help to define the probe position onto a boundary as expected
* by the sdis_solve_probe_boundary function.
@@ -912,6 +945,11 @@ sdis_estimator_get_total_flux
struct sdis_mc* flux);
SDIS_API res_T
+sdis_estimator_get_power
+ (const struct sdis_estimator* estimator,
+ struct sdis_mc* power);
+
+SDIS_API res_T
sdis_estimator_get_paths_count
(const struct sdis_estimator* estimator,
size_t* npaths);
@@ -953,6 +991,17 @@ sdis_green_function_solve
const double time_range[2], /* Observation time */
struct sdis_estimator** estimator);
+SDIS_API res_T
+sdis_green_function_write
+ (struct sdis_green_function* green,
+ FILE* stream);
+
+SDIS_API res_T
+sdis_green_function_create_from_stream
+ (struct sdis_scene* scn, /* Scene from which the green was evaluated */
+ FILE* stream, /* Stream into which the green was serialized */
+ struct sdis_green_function** green);
+
/* Retrieve the number of valid paths used to estimate the green function. It
* is actually equal to the number of successful realisations. */
SDIS_API res_T
@@ -1083,6 +1132,15 @@ sdis_solve_medium
const struct sdis_solve_medium_args* args,
struct sdis_estimator** estimator);
+/* P = SUM(volumic_power(x)) / Nrealisations * Volume
+ * power (in Watt) = time_range[0] == time_range[1]
+ * ? P : P / (time_range[1] - time_range[0]) */
+SDIS_API res_T
+sdis_compute_power
+ (struct sdis_scene* scn,
+ const struct sdis_compute_power_args* args,
+ struct sdis_estimator** estimator);
+
/*******************************************************************************
* Green solvers.
*
diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h
@@ -27,7 +27,7 @@
struct ssp_rng;
struct ssp_rng_proxy;
-struct name { FITEM; };
+struct name { FITEM; void* mem; };
#define FITEM_TYPE name
#include <rsys/free_list.h>
diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c
@@ -98,7 +98,10 @@ res_T
sdis_estimator_get_temperature
(const struct sdis_estimator* estimator, struct sdis_mc* mc)
{
- if(!estimator || !mc) return RES_BAD_ARG;
+ if(!estimator || !mc
+ || ( estimator->type != SDIS_ESTIMATOR_TEMPERATURE
+ && estimator->type != SDIS_ESTIMATOR_FLUX))
+ return RES_BAD_ARG;
SETUP_MC(mc, &estimator->temperature);
return RES_OK;
}
@@ -116,9 +119,8 @@ res_T
sdis_estimator_get_convective_flux
(const struct sdis_estimator* estimator, struct sdis_mc* flux)
{
- if(!estimator || !flux ||estimator->type != SDIS_ESTIMATOR_FLUX)
+ if(!estimator || !flux || estimator->type != SDIS_ESTIMATOR_FLUX)
return RES_BAD_ARG;
- ASSERT(estimator->fluxes);
SETUP_MC(flux, &estimator->fluxes[FLUX_CONVECTIVE]);
return RES_OK;
}
@@ -129,7 +131,6 @@ sdis_estimator_get_radiative_flux
{
if(!estimator || !flux || estimator->type != SDIS_ESTIMATOR_FLUX)
return RES_BAD_ARG;
- ASSERT(estimator->fluxes);
SETUP_MC(flux, &estimator->fluxes[FLUX_RADIATIVE]);
return RES_OK;
}
@@ -140,11 +141,27 @@ sdis_estimator_get_total_flux
{
if(!estimator || !flux || estimator->type != SDIS_ESTIMATOR_FLUX)
return RES_BAD_ARG;
- ASSERT(estimator->fluxes);
SETUP_MC(flux, &estimator->fluxes[FLUX_TOTAL]);
return RES_OK;
}
+res_T
+sdis_estimator_get_power
+ (const struct sdis_estimator* estimator, struct sdis_mc* power)
+{
+ if(!estimator || !power || estimator->type != SDIS_ESTIMATOR_POWER)
+ return RES_BAD_ARG;
+ SETUP_MC(power, &estimator->power.power);
+ power->E *= estimator->power.spread;
+
+ if(estimator->power.time_range[0]
+ != estimator->power.time_range[1]) {
+ power->E /= /* From Joule to Watt */
+ (estimator->power.time_range[1]-estimator->power.time_range[0]);
+ }
+ return RES_OK;
+}
+
#undef SETUP_MC
res_T
diff --git a/src/sdis_estimator_c.h b/src/sdis_estimator_c.h
@@ -39,6 +39,13 @@ struct sdis_estimator {
struct accum temperature;
struct accum realisation_time;
struct accum fluxes[FLUX_NAMES_COUNT__];
+
+ struct {
+ struct accum power;
+ double spread;
+ double time_range[2];
+ } power;
+
size_t nrealisations; /* #successes */
size_t nfailures;
@@ -98,6 +105,23 @@ estimator_setup_temperature
}
static INLINE void
+estimator_setup_power
+ (struct sdis_estimator* estim,
+ const double sum,
+ const double sum2,
+ const double spread,
+ const double time_range[2])
+{
+ ASSERT(estim && estim->nrealisations && time_range);
+ estim->power.power.sum = sum;
+ estim->power.power.sum2 = sum2;
+ estim->power.power.count = estim->nrealisations;
+ estim->power.spread = spread;
+ estim->power.time_range[0] = time_range[0];
+ estim->power.time_range[1] = time_range[1];
+}
+
+static INLINE void
estimator_setup_realisation_time
(struct sdis_estimator* estim,
const double sum,
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -16,10 +16,11 @@
#include "sdis_device_c.h"
#include "sdis_estimator_c.h"
#include "sdis_green.h"
+#include "sdis_interface_c.h"
#include "sdis_log.h"
#include "sdis_medium_c.h"
#include "sdis_misc.h"
-#include "sdis_interface_c.h"
+#include "sdis_scene_c.h"
#include <star/ssp.h>
@@ -98,6 +99,7 @@ green_path_init(struct mem_allocator* allocator, struct green_path* path)
darray_flux_term_init(allocator, &path->flux_terms);
darray_power_term_init(allocator, &path->power_terms);
path->limit.vertex = SDIS_RWALK_VERTEX_NULL;
+ path->limit.fragment = SDIS_INTERFACE_FRAGMENT_NULL;
path->limit_id = UINT_MAX;
path->limit_type = SDIS_POINT_NONE;
path->ilast_medium = UINT16_MAX;
@@ -164,6 +166,96 @@ green_path_copy_and_release(struct green_path* dst, struct green_path* src)
return RES_OK;
}
+static res_T
+green_path_write(const struct green_path* path, FILE* stream)
+{
+ size_t sz = 0;
+ res_T res = RES_OK;
+ ASSERT(path && stream);
+
+ #define WRITE(Var, N) { \
+ if(fwrite((Var), sizeof(*(Var)), (N), stream) != (N)) { \
+ res = RES_IO_ERR; \
+ goto error; \
+ } \
+ } (void)0
+
+ /* Write the list of flux terms */
+ sz = darray_flux_term_size_get(&path->flux_terms);
+ WRITE(&sz, 1);
+ WRITE(darray_flux_term_cdata_get(&path->flux_terms), sz);
+
+ /* Write the list of power terms */
+ sz = darray_power_term_size_get(&path->power_terms);
+ WRITE(&sz, 1);
+ WRITE(darray_power_term_cdata_get(&path->power_terms), sz);
+
+ /* Write the limit point */
+ WRITE(&path->limit, 1);
+ WRITE(&path->limit_id, 1);
+ WRITE(&path->limit_type, 1);
+
+ /* Write miscellaneous data */
+ WRITE(&path->ilast_medium, 1);
+ WRITE(&path->ilast_interf, 1);
+
+ #undef WRITE
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+green_path_read(struct green_path* path, FILE* stream)
+{
+ size_t sz = 0;
+ res_T res = RES_OK;
+ ASSERT(path && stream);
+
+ #define READ(Var, N) { \
+ if(fread((Var), sizeof(*(Var)), (N), stream) != (N)) { \
+ if(feof(stream)) { \
+ res = RES_BAD_ARG; \
+ } else if(ferror(stream)) { \
+ res = RES_IO_ERR; \
+ } else { \
+ res = RES_UNKNOWN_ERR; \
+ } \
+ goto error; \
+ } \
+ } (void)0
+
+ /* Read the list of flux terms */
+ READ(&sz, 1);
+ res = darray_flux_term_resize(&path->flux_terms, sz);
+ if(res != RES_OK) goto error;
+ READ(darray_flux_term_data_get(&path->flux_terms), sz);
+
+ /* Read the list of power tems */
+ READ(&sz, 1);
+ res = darray_power_term_resize(&path->power_terms, sz);
+ if(res != RES_OK) goto error;
+ READ(darray_power_term_data_get(&path->power_terms), sz);
+
+ /* Read the limit point */
+ READ(&path->limit, 1);
+ READ(&path->limit_id, 1);
+ READ(&path->limit_type, 1);
+
+ /* Read the miscellaneous data */
+ READ(&path->ilast_medium, 1);
+ READ(&path->ilast_interf, 1);
+
+ #undef READ
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/* Generate the dynamic array of green paths */
#define DARRAY_NAME green_path
#define DARRAY_DATA struct green_path
@@ -199,12 +291,44 @@ struct sdis_green_function {
FILE* rng_state;
ref_T ref;
- struct sdis_device* dev;
+ struct sdis_scene* scn;
};
/*******************************************************************************
* Helper functions
******************************************************************************/
+enum rng_type {
+ RNG_KISS,
+ RNG_MT_19937_64,
+ RNG_RANLUX48,
+ RNG_THREEFRY,
+ RNG_UNKNOWN
+};
+
+static INLINE enum rng_type
+get_rng_type_enum(const struct ssp_rng_type* type)
+{
+ ASSERT(type);
+ if(ssp_rng_type_eq(type, &ssp_rng_kiss)) return RNG_KISS;
+ if(ssp_rng_type_eq(type, &ssp_rng_mt19937_64)) return RNG_MT_19937_64;
+ if(ssp_rng_type_eq(type, &ssp_rng_ranlux48)) return RNG_RANLUX48;
+ if(ssp_rng_type_eq(type, &ssp_rng_threefry)) return RNG_THREEFRY;
+ return RNG_UNKNOWN;
+}
+
+static INLINE void
+get_rng_type_ssp(const enum rng_type type, struct ssp_rng_type* ssp_type)
+{
+ switch(type) {
+ case RNG_KISS: *ssp_type = ssp_rng_kiss; break;
+ case RNG_MT_19937_64: *ssp_type = ssp_rng_mt19937_64; break;
+ case RNG_RANLUX48: *ssp_type = ssp_rng_ranlux48; break;
+ case RNG_THREEFRY: *ssp_type = ssp_rng_threefry; break;
+ case RNG_UNKNOWN: memset(ssp_type, 0, sizeof(*ssp_type)); break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+}
+
static res_T
ensure_medium_registration
(struct sdis_green_function* green,
@@ -338,7 +462,7 @@ green_function_solve_path
if(time_curr <= 0
|| (path->limit_type == SDIS_VERTEX && time_curr <= medium_get_t0(medium))) {
- log_err(green->dev,
+ log_err(green->scn->dev,
"%s: invalid observation time \"%g\": the initial condition is reached "
"while instationary system are not supported by the green function.\n",
FUNC_NAME, time);
@@ -370,6 +494,280 @@ error:
goto exit;
}
+static res_T
+write_media(struct sdis_green_function* green, FILE* stream)
+{
+ struct htable_medium_iterator it, it_end;
+ size_t nmedia = 0;
+ res_T res = RES_OK;
+ ASSERT(green && stream);
+
+ #define WRITE(Var) { \
+ if(fwrite((Var), sizeof(*(Var)), 1, stream) != 1) { \
+ res = RES_IO_ERR; \
+ goto error; \
+ } \
+ } (void)0
+
+ nmedia = htable_medium_size_get(&green->media);
+ WRITE(&nmedia);
+
+ htable_medium_begin(&green->media, &it);
+ htable_medium_end(&green->media, &it_end);
+ while(!htable_medium_iterator_eq(&it, &it_end)) {
+ const struct sdis_medium* mdm = *htable_medium_iterator_data_get(&it);
+ htable_medium_iterator_next(&it);
+ WRITE(&mdm->id);
+ WRITE(&mdm->type);
+ }
+
+ #undef WRITE
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+read_media(struct sdis_green_function* green, FILE* stream)
+{
+ size_t nmedia = 0;
+ size_t imedium = 0;
+ res_T res = RES_OK;
+ ASSERT(green && stream);
+
+ #define READ(Var) { \
+ if(fread((Var), sizeof(*(Var)), 1, stream) != 1) { \
+ if(feof(stream)) { \
+ res = RES_BAD_ARG; \
+ } else if(ferror(stream)) { \
+ res = RES_IO_ERR; \
+ } else { \
+ res = RES_UNKNOWN_ERR; \
+ } \
+ goto error; \
+ } \
+ } (void)0
+
+ READ(&nmedia);
+ FOR_EACH(imedium, 0, nmedia) {
+ struct name* name = NULL;
+ struct sdis_medium* mdm = NULL;
+ struct fid id;
+ enum sdis_medium_type mdm_type;
+
+ READ(&id);
+ READ(&mdm_type);
+
+ name = flist_name_get(&green->scn->dev->media_names, id);
+ if(!name) {
+ log_err(green->scn->dev, "%s: a Stardis medium is missing.\n",
+ FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ mdm = name->mem;
+
+ if(mdm_type != mdm->type) {
+ log_err(green->scn->dev, "%s: inconsistency between the a Stardis medium "
+ "and its serialised data.\n", FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = ensure_medium_registration(green, mdm);
+ if(res != RES_OK) goto error;
+ }
+
+ #undef READ
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+write_interfaces(struct sdis_green_function* green, FILE* stream)
+{
+ struct htable_interf_iterator it, it_end;
+ size_t ninterfaces = 0;
+ res_T res = RES_OK;
+ ASSERT(green && stream);
+
+ #define WRITE(Var) { \
+ if(fwrite((Var), sizeof(*(Var)), 1, stream) != 1) { \
+ res = RES_IO_ERR; \
+ goto error; \
+ } \
+ } (void)0
+ ninterfaces = htable_interf_size_get(&green->interfaces);
+ WRITE(&ninterfaces);
+
+ htable_interf_begin(&green->interfaces, &it);
+ htable_interf_end(&green->interfaces, &it_end);
+ while(!htable_interf_iterator_eq(&it, &it_end)) {
+ const struct sdis_interface* interf = *htable_interf_iterator_data_get(&it);
+ htable_interf_iterator_next(&it);
+ WRITE(&interf->id);
+ WRITE(&interf->medium_front->id);
+ WRITE(&interf->medium_front->type);
+ WRITE(&interf->medium_back->id);
+ WRITE(&interf->medium_back->type);
+ }
+ #undef WRITE
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+read_interfaces(struct sdis_green_function* green, FILE* stream)
+{
+ size_t ninterfs = 0;
+ size_t iinterf = 0;
+ res_T res = RES_OK;
+ ASSERT(green && stream);
+
+ #define READ(Var) { \
+ if(fread((Var), sizeof(*(Var)), 1, stream) != 1) { \
+ if(feof(stream)) { \
+ res = RES_BAD_ARG; \
+ } else if(ferror(stream)) { \
+ res = RES_IO_ERR; \
+ } else { \
+ res = RES_UNKNOWN_ERR; \
+ } \
+ goto error; \
+ } \
+ } (void)0
+
+ READ(&ninterfs);
+ FOR_EACH(iinterf, 0, ninterfs) {
+ struct name* name = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* mdm_front = NULL;
+ struct sdis_medium* mdm_back = NULL;
+ struct fid id;
+ struct fid mdm_front_id;
+ struct fid mdm_back_id;
+ enum sdis_medium_type mdm_front_type;
+ enum sdis_medium_type mdm_back_type;
+
+ READ(&id);
+ READ(&mdm_front_id);
+ READ(&mdm_front_type);
+ READ(&mdm_back_id);
+ READ(&mdm_back_type);
+
+ name = flist_name_get(&green->scn->dev->interfaces_names, id);
+ if(!name) {
+ log_err(green->scn->dev, "%s: a Stardis interface is missing.\n",
+ FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ interf = name->mem;
+ mdm_front = flist_name_get(&green->scn->dev->media_names, mdm_front_id)->mem;
+ mdm_back = flist_name_get(&green->scn->dev->media_names, mdm_back_id)->mem;
+
+ if(mdm_front != interf->medium_front
+ || mdm_back != interf->medium_back
+ || mdm_front_type != interf->medium_front->type
+ || mdm_back_type != interf->medium_back->type) {
+ log_err(green->scn->dev, "%s: inconsistency between the a Stardis interface "
+ "and its serialised data.\n", FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = ensure_interface_registration(green, interf);
+ if(res != RES_OK) goto error;
+ }
+
+ #undef READ
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+write_paths_list(struct sdis_green_function* green, FILE* stream)
+{
+ size_t npaths = 0;
+ size_t ipath = 0;
+ res_T res = RES_OK;
+ ASSERT(green && stream);
+
+ #define WRITE(Var) { \
+ if(fwrite((Var), sizeof(*(Var)), 1, stream) != 1) { \
+ res = RES_IO_ERR; \
+ goto error; \
+ } \
+ } (void)0
+ npaths = darray_green_path_size_get(&green->paths);
+ WRITE(&npaths);
+ FOR_EACH(ipath, 0, npaths) {
+ const struct green_path* path = NULL;
+ path = darray_green_path_cdata_get(&green->paths) + ipath;
+
+ res = green_path_write(path, stream);
+ if(res != RES_OK) goto error;
+ }
+ #undef WRITE
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+read_paths_list(struct sdis_green_function* green, FILE* stream)
+{
+ size_t npaths = 0;
+ size_t ipath = 0;
+ res_T res = RES_OK;
+
+ #define READ(Var) { \
+ if(fread((Var), sizeof(*(Var)), 1, stream) != 1) { \
+ if(feof(stream)) { \
+ res = RES_BAD_ARG; \
+ } else if(ferror(stream)) { \
+ res = RES_IO_ERR; \
+ } else { \
+ res = RES_UNKNOWN_ERR; \
+ } \
+ goto error; \
+ } \
+ } (void)0
+
+ READ(&npaths);
+ res = darray_green_path_resize(&green->paths, npaths);
+ if(res != RES_OK) goto error;
+
+ FOR_EACH(ipath, 0, npaths) {
+ struct green_path* path = NULL;
+ path = darray_green_path_data_get(&green->paths) + ipath;
+
+ res = green_path_read(path, stream);
+ if(res != RES_OK) goto error;
+ }
+ #undef READ
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
static void
green_function_clear(struct sdis_green_function* green)
{
@@ -406,18 +804,18 @@ green_function_clear(struct sdis_green_function* green)
static void
green_function_release(ref_T* ref)
{
- struct sdis_device* dev;
+ struct sdis_scene* scn;
struct sdis_green_function* green;
ASSERT(ref);
green = CONTAINER_OF(ref, struct sdis_green_function, ref);
- dev = green->dev;
+ scn = green->scn;
green_function_clear(green);
htable_medium_release(&green->media);
htable_interf_release(&green->interfaces);
darray_green_path_release(&green->paths);
if(green->rng_state) fclose(green->rng_state);
- MEM_RM(dev->allocator, green);
- SDIS(device_ref_put(dev));
+ MEM_RM(scn->dev->allocator, green);
+ SDIS(scene_ref_put(scn));
}
/*******************************************************************************
@@ -460,7 +858,7 @@ sdis_green_function_solve
goto error;
}
- res = ssp_rng_create(green->dev->allocator, &green->rng_type, &rng);
+ res = ssp_rng_create(green->scn->dev->allocator, &green->rng_type, &rng);
if(res != RES_OK) goto error;
/* Avoid correlation by defining the RNG state from the final state of the
@@ -472,7 +870,7 @@ sdis_green_function_solve
npaths = darray_green_path_size_get(&green->paths);
/* Create the estimator */
- res = estimator_create(green->dev, SDIS_ESTIMATOR_TEMPERATURE, &estimator);
+ res = estimator_create(green->scn->dev, SDIS_ESTIMATOR_TEMPERATURE, &estimator);
if(res != RES_OK) goto error;
/* Solve the green function */
@@ -508,6 +906,161 @@ error:
}
res_T
+sdis_green_function_write(struct sdis_green_function* green, FILE* stream)
+{
+ struct ssp_rng* rng = NULL;
+ enum rng_type rng_type = RNG_UNKNOWN;
+ hash256_T hash;
+ res_T res = RES_OK;
+
+ if(!green || !stream) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ #define WRITE(Var, Nb) { \
+ if(fwrite((Var), sizeof(*(Var)), (Nb), stream) != (Nb)) { \
+ res = RES_IO_ERR; \
+ goto error; \
+ } \
+ } (void)0
+
+ rng_type = get_rng_type_enum(&green->rng_type);
+ if(rng_type == RNG_UNKNOWN) {
+ log_err(green->scn->dev,
+ "%s: could not function a green function with an unknown RNG type.\n",
+ FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ WRITE(&SDIS_GREEN_FUNCTION_VERSION, 1);
+
+ res = scene_compute_hash(green->scn, hash);
+ if(res != RES_OK) goto error;
+ WRITE(hash, sizeof(hash256_T));
+
+ res = write_media(green, stream);
+ if(res != RES_OK) goto error;
+ res = write_interfaces(green, stream);
+ if(res != RES_OK) goto error;
+ res = write_paths_list(green, stream);
+ if(res != RES_OK) goto error;
+
+ WRITE(&green->npaths_valid, 1);
+ WRITE(&green->npaths_invalid, 1);
+ WRITE(&green->realisation_time, 1);
+ WRITE(&rng_type, 1);
+ #undef WRITE
+
+ /* Create a temporary RNG used to serialise the RNG state */
+ res = ssp_rng_create(green->scn->dev->allocator, &green->rng_type, &rng);
+ if(res != RES_OK) goto error;
+ rewind(green->rng_state);
+ res = ssp_rng_read(rng, green->rng_state);
+ if(res != RES_OK) goto error;
+ res = ssp_rng_write(rng, stream);
+ if(res != RES_OK) goto error;
+
+exit:
+ if(rng) SSP(rng_ref_put(rng));
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+sdis_green_function_create_from_stream
+ (struct sdis_scene* scn,
+ FILE* stream,
+ struct sdis_green_function** out_green)
+{
+ hash256_T hash0, hash1;
+ struct sdis_green_function* green = NULL;
+ struct ssp_rng* rng = NULL;
+ enum rng_type rng_type = RNG_UNKNOWN;
+ int version = 0;
+ res_T res = RES_OK;
+
+ if(!scn || !stream || !out_green) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = green_function_create(scn, &green);
+ if(res != RES_OK) goto error;
+
+ #define READ(Var, Nb) { \
+ if(fread((Var), sizeof(*(Var)), (Nb), stream) != (Nb)) { \
+ if(feof(stream)) { \
+ res = RES_BAD_ARG; \
+ } else if(ferror(stream)) { \
+ res = RES_IO_ERR; \
+ } else { \
+ res = RES_UNKNOWN_ERR; \
+ } \
+ goto error; \
+ } \
+ } (void)0
+
+ READ(&version, 1);
+ if(version != SDIS_GREEN_FUNCTION_VERSION) {
+ log_err(green->scn->dev,
+ "%s: unexpected green function version %d. Expecting a green function "
+ "in version %d.\n",
+ FUNC_NAME, version, SDIS_GREEN_FUNCTION_VERSION);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = scene_compute_hash(green->scn, hash0);
+ if(res != RES_OK) goto error;
+
+ READ(hash1, sizeof(hash256_T));
+ if(!hash256_eq(hash0, hash1)) {
+ log_err(green->scn->dev,
+ "%s: the submitted scene does not match scene used to estimate the green "
+ "function.\n", FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = read_media(green, stream);
+ if(res != RES_OK) goto error;
+ res = read_interfaces(green, stream);
+ if(res != RES_OK) goto error;
+ res = read_paths_list(green, stream);
+ if(res != RES_OK) goto error;
+
+ READ(&green->npaths_valid, 1);
+ READ(&green->npaths_invalid, 1);
+ READ(&green->realisation_time, 1);
+ READ(&rng_type, 1);
+ #undef READ
+
+ get_rng_type_ssp(rng_type, &green->rng_type);
+
+ /* Create a temporary RNG used to deserialise the RNG state */
+ res = ssp_rng_create(green->scn->dev->allocator, &green->rng_type, &rng);
+ if(res != RES_OK) goto error;
+ res = ssp_rng_read(rng, stream);
+ if(res != RES_OK) goto error;
+ res = ssp_rng_write(rng, green->rng_state);
+ if(res != RES_OK) goto error;
+
+exit:
+ if(rng) SSP(rng_ref_put(rng));
+ if(out_green) *out_green = green;
+ return res;
+error:
+ if(green) {
+ SDIS(green_function_ref_put(green));
+ green = NULL;
+ }
+ goto exit;
+}
+
+res_T
sdis_green_function_get_paths_count
(const struct sdis_green_function* green, size_t* npaths)
{
@@ -733,23 +1286,23 @@ error:
******************************************************************************/
res_T
green_function_create
- (struct sdis_device* dev, struct sdis_green_function** out_green)
+ (struct sdis_scene* scn, struct sdis_green_function** out_green)
{
struct sdis_green_function* green = NULL;
res_T res = RES_OK;
- ASSERT(dev && out_green);
+ ASSERT(scn && out_green);
- green = MEM_CALLOC(dev->allocator, 1, sizeof(*green));
+ green = MEM_CALLOC(scn->dev->allocator, 1, sizeof(*green));
if(!green) {
res = RES_MEM_ERR;
goto error;
}
ref_init(&green->ref);
- SDIS(device_ref_get(dev));
- green->dev = dev;
- htable_medium_init(dev->allocator, &green->media);
- htable_interf_init(dev->allocator, &green->interfaces);
- darray_green_path_init(dev->allocator, &green->paths);
+ SDIS(scene_ref_get(scn));
+ green->scn = scn;
+ htable_medium_init(scn->dev->allocator, &green->media);
+ htable_interf_init(scn->dev->allocator, &green->interfaces);
+ darray_green_path_init(scn->dev->allocator, &green->paths);
green->npaths_valid = SIZE_MAX;
green->npaths_invalid = SIZE_MAX;
diff --git a/src/sdis_green.h b/src/sdis_green.h
@@ -18,6 +18,11 @@
#include <rsys/rsys.h>
+/* Current version the green function data structure. One should increment it
+ * and perform a version management onto serialized data when the gren function
+ * data structure is updated. */
+static const int SDIS_GREEN_FUNCTION_VERSION = 0;
+
/* Forward declaration */
struct accum;
struct sdis_green_function;
@@ -34,7 +39,7 @@ static const struct green_path_handle GREEN_PATH_HANDLE_NULL =
extern LOCAL_SYM res_T
green_function_create
- (struct sdis_device* dev,
+ (struct sdis_scene* scn,
struct sdis_green_function** green);
/* Merge `src' into `dst' an clear `src' */
diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h
@@ -600,6 +600,11 @@ XD(solid_solid_boundary_path)
}
}
+ /* Time rewind */
+ res = XD(time_rewind)(mdm, rng, reinject_dst, fp_to_meter, ctx, rwalk, T);
+ if(res != RES_OK) goto error;
+ if(T->done) goto exit; /* Limit condition was reached */
+
/* Perform reinjection. */
XD(move_pos)(rwalk->vtx.P, dir, (float)reinject_dst);
if(hit->distance == reinject_dst) {
@@ -764,6 +769,11 @@ XD(solid_fluid_boundary_path)
}
}
+ /* Time rewind */
+ res = XD(time_rewind)(solid, rng, reinject_dst, fp_to_meter, ctx, rwalk, T);
+ if(res != RES_OK) goto error;
+ if(T->done) goto exit; /* Limit condition was reached */
+
/* Perform solid reinjection */
XD(move_pos)(rwalk->vtx.P, dir0, reinject_dst);
if(hit.distance == reinject_dst) {
@@ -898,6 +908,11 @@ XD(solid_boundary_with_flux_path)
}
}
+ /* Time rewind */
+ res = XD(time_rewind)(mdm, rng, reinject_dst, fp_to_meter, ctx, rwalk, T);
+ if(res != RES_OK) goto error;
+ if(T->done) goto exit; /* Limit condition was reached */
+
/* Reinject. If the reinjection move the point too close of a boundary,
* assume that the zone is isotherm and move to the boundary. */
XD(move_pos)(rwalk->vtx.P, dir0, reinject_dst);
@@ -911,7 +926,6 @@ XD(solid_boundary_with_flux_path)
rwalk->mdm = mdm;
rwalk->hit = SXD_HIT_NULL;
rwalk->hit_side = SDIS_SIDE_NULL__;
-
}
/* Register the new vertex against the heat path */
@@ -975,7 +989,7 @@ XD(boundary_path)
/* Check if the boundary flux is known. Note that currently, only solid media
* can have a flux as limit condition */
mdm = interface_get_medium(interf, frag.side);
- if(sdis_medium_get_type(mdm) == SDIS_SOLID ) {
+ if(sdis_medium_get_type(mdm) == SDIS_SOLID) {
const double phi = interface_side_get_flux(interf, &frag);
if(phi != SDIS_FLUX_NONE) {
res = XD(solid_boundary_with_flux_path)
diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h
@@ -235,8 +235,6 @@ XD(conductive_path)
do { /* Solid random walk */
struct sXd(hit) hit0, hit1;
double lambda; /* Thermal conductivity */
- double rho; /* Volumic mass */
- double cp; /* Calorific capacity */
double tmp;
double power_factor = 0;
double power;
@@ -244,7 +242,8 @@ XD(conductive_path)
float dir0[DIM], dir1[DIM];
float org[DIM];
- /* Check the limit condition */
+ /* Check the limit condition
+ * REVIEW Rfo: This can be a bug if the random walk comes from a boundary */
tmp = solid_get_temperature(mdm, &rwalk->vtx);
if(tmp >= 0) {
T->value += tmp;
@@ -266,8 +265,6 @@ XD(conductive_path)
/* Fetch solid properties */
delta_solid = (float)solid_get_delta(mdm, &rwalk->vtx);
lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
- rho = solid_get_volumic_mass(mdm, &rwalk->vtx);
- cp = solid_get_calorific_capacity(mdm, &rwalk->vtx);
power = solid_get_volumic_power(mdm, &rwalk->vtx);
if(ctx->green_path && power_ref != power) {
@@ -345,39 +342,12 @@ XD(conductive_path)
green_power_factor += power_factor;
}
- /* Sample the time */
- if(!IS_INF(rwalk->vtx.time)) {
- double tau, mu, t0;
- mu = (2*DIM*lambda) / (rho*cp*delta*fp_to_meter*delta*fp_to_meter);
- tau = ssp_ran_exp(rng, mu);
- t0 = ctx->green_path ? -INF : solid_get_t0(rwalk->mdm);
- rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, t0);
- if(rwalk->vtx.time == t0) {
- /* Check the initial condition */
- tmp = solid_get_temperature(mdm, &rwalk->vtx);
- if(tmp >= 0) {
- T->value += tmp;
- T->done = 1;
-
- if(ctx->heat_path) {
- struct sdis_heat_vertex* vtx;
- vtx = heat_path_get_last_vertex(ctx->heat_path);
- vtx->time = rwalk->vtx.time;
- vtx->weight = T->value;
- }
- break;
- }
- /* The initial condition should have been reached */
- log_err(scn->dev,
- "%s: undefined initial condition. "
- "The time is %g but the temperature remains unknown.\n",
- FUNC_NAME, t0);
- res = RES_BAD_OP;
- goto error;
- }
- }
+ /* Rewind the time */
+ res = XD(time_rewind)(rwalk->mdm, rng, delta, fp_to_meter, ctx, rwalk, T);
+ if(res != RES_OK) goto error;
+ if(T->done) break; /* Limit condition was reached */
- /* Define if the random walk hits something along dir0 */
+ /* Define if the random walk hits something along dir0 */
if(hit0.distance > delta) {
rwalk->hit = SXD_HIT_NULL;
rwalk->hit_side = SDIS_SIDE_NULL__;
diff --git a/src/sdis_interface.c b/src/sdis_interface.c
@@ -153,6 +153,7 @@ sdis_interface_create
interf->dev = dev;
interf->shader = *shader;
interf->id = flist_name_add(&dev->interfaces_names);
+ flist_name_get(&dev->interfaces_names, interf->id)->mem = interf;
if(data) {
SDIS(data_ref_get(data));
diff --git a/src/sdis_medium.c b/src/sdis_medium.c
@@ -68,6 +68,7 @@ medium_create
medium->dev = dev;
medium->type = type;
medium->id = flist_name_add(&dev->media_names);
+ flist_name_get(&dev->media_names, medium->id)->mem = medium;
exit:
if(out_medium) *out_medium = medium;
diff --git a/src/sdis_medium_c.h b/src/sdis_medium_c.h
@@ -18,7 +18,9 @@
#include "sdis.h"
+#include <rsys/free_list.h>
#include <rsys/math.h>
+#include <rsys/ref_count.h>
struct sdis_medium {
enum sdis_medium_type type;
diff --git a/src/sdis_misc.c b/src/sdis_misc.c
@@ -0,0 +1,23 @@
+/* Copyright (C) 2016-2020 |Meso|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"
+
+/* Generate the generic functions */
+#define SDIS_XD_DIMENSION 2
+#include "sdis_misc_Xd.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_misc_Xd.h"
+
diff --git a/src/sdis_misc.h b/src/sdis_misc.h
@@ -134,4 +134,24 @@ register_heat_vertex
return heat_path_add_vertex(path, &heat_vtx);
}
+extern LOCAL_SYM res_T
+time_rewind_2d
+ (const struct sdis_medium* mdm, /* Medium into which the time is rewinded */
+ struct ssp_rng* rng,
+ const double delta,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct rwalk_2d* rwalk,
+ struct temperature_2d* T);
+
+extern LOCAL_SYM res_T
+time_rewind_3d
+ (const struct sdis_medium* mdm, /* Medium into which the time is rewinded */
+ struct ssp_rng* rng,
+ const double delta,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct rwalk_3d* rwalk,
+ struct temperature_3d* T);
+
#endif /* SDIS_MISC_H */
diff --git a/src/sdis_misc_Xd.h b/src/sdis_misc_Xd.h
@@ -0,0 +1,91 @@
+/* Copyright (C) 2016-2020 |Meso|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_heat_path.h"
+#include "sdis_log.h"
+#include "sdis_medium_c.h"
+#include "sdis_misc.h"
+
+#include <star/ssp.h>
+
+#include "sdis_Xd_begin.h"
+
+res_T
+XD(time_rewind)
+ (const struct sdis_medium* mdm,
+ struct ssp_rng* rng,
+ const double delta,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct XD(temperature)* T)
+{
+ const double delta_in_meter = delta * fp_to_meter;
+ double temperature;
+ double lambda, rho, cp;
+ double tau, mu, t0;
+ res_T res = RES_OK;
+ ASSERT(mdm && rng && delta && fp_to_meter && ctx && rwalk);
+ ASSERT(sdis_medium_get_type(mdm) == SDIS_SOLID);
+ ASSERT(T->done == 0);
+
+ if(IS_INF(rwalk->vtx.time)) goto exit;
+
+ /* Fetch phyisical properties */
+ lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
+ rho = solid_get_volumic_mass(mdm, &rwalk->vtx);
+ cp = solid_get_calorific_capacity(mdm, &rwalk->vtx);
+
+ /* Fetch the limit time */
+ t0 = ctx->green_path ? -INF : solid_get_t0(mdm);
+
+ /* Sample the time to reroll */
+ mu = (2*DIM*lambda)/(rho*cp*delta_in_meter*delta_in_meter);
+ tau = ssp_ran_exp(rng, mu);
+
+ /* Time rewind */
+ rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, t0);
+
+ /* The path does not reach the limit condition */
+ if(rwalk->vtx.time > t0) goto exit;
+
+ /* Fetch initial temperature */
+ temperature = solid_get_temperature(mdm, &rwalk->vtx);
+ if(temperature < 0) {
+ log_err(mdm->dev, "%s: the path reaches the limit condition by the "
+ "temperature remains unknown.\n", FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Update temperature */
+ T->value += temperature;
+ T->done = 1;
+
+ if(ctx->heat_path) {
+ /* Update the registered vertex data */
+ struct sdis_heat_vertex* vtx;
+ vtx = heat_path_get_last_vertex(ctx->heat_path);
+ vtx->time = rwalk->vtx.time;
+ vtx->weight = T->value;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#include "sdis_Xd_end.h"
diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h
@@ -144,4 +144,3 @@ ray_realisation_3d
double* weight);
#endif /* SDIS_REALISATION_H */
-
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -13,6 +13,10 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
+#define _POSIX_C_SOURCE 200809L /* mmap support */
+#define _DEFAULT_SOURCE 1 /* MAP_POPULATE support */
+#define _BSD_SOURCE 1 /* MAP_POPULATE for glibc < 2.19 */
+
#include "sdis_scene_Xd.h"
/* Generate the Generic functions of the scene */
@@ -22,9 +26,12 @@
#include "sdis_scene_Xd.h"
#include "sdis.h"
+#include "sdis_interface_c.h"
#include "sdis_scene_c.h"
+#include <float.h>
#include <limits.h>
+#include <sys/mman.h>
/*******************************************************************************
* Helper function
@@ -182,6 +189,22 @@ sdis_scene_get_aabb
}
res_T
+sdis_scene_find_closest_point
+ (const struct sdis_scene* scn,
+ const double pos[3],
+ const double radius,
+ size_t* iprim,
+ double uv[2])
+{
+ if(!scn) return RES_BAD_ARG;
+ if(scene_is_2d(scn)) {
+ return scene_find_closest_point_2d(scn, pos, radius, iprim, uv);
+ } else {
+ return scene_find_closest_point_3d(scn, pos, radius, iprim, uv);
+ }
+}
+
+res_T
sdis_scene_get_boundary_position
(const struct sdis_scene* scn,
const size_t iprim,
@@ -333,7 +356,7 @@ sdis_scene_get_medium_spread
}
}
*out_spread = spread;
-
+
exit:
return res;
error:
@@ -373,3 +396,93 @@ scene_get_medium_in_closed_boundaries
: scene_get_medium_in_closed_boundaries_3d(scn, pos, out_medium);
}
+res_T
+scene_compute_hash(const struct sdis_scene* scn, hash256_T hash)
+{
+ void* data = NULL;
+ FILE* stream = NULL;
+ size_t iprim, nprims;
+ size_t len;
+ res_T res = RES_OK;
+ ASSERT(scn && hash);
+
+ stream = tmpfile();
+ if(!stream) {
+ res = RES_IO_ERR;
+ goto error;
+ }
+
+ #define WRITE(Var, Nb) \
+ if(fwrite((Var), sizeof(*(Var)), Nb, stream) != (Nb)) { \
+ res = RES_IO_ERR; \
+ goto error; \
+ } (void)0
+ if(scene_is_2d(scn)) {
+ S2D(scene_view_primitives_count(scn->s2d_view, &nprims));
+ } else {
+ S3D(scene_view_primitives_count(scn->s3d_view, &nprims));
+ }
+ FOR_EACH(iprim, 0, nprims) {
+ struct sdis_interface* interf = NULL;
+ size_t ivert;
+
+ if(scene_is_2d(scn)) {
+ struct s2d_primitive prim;
+ S2D(scene_view_get_primitive(scn->s2d_view, (unsigned)iprim, &prim));
+ FOR_EACH(ivert, 0, 2) {
+ struct s2d_attrib attr;
+ S2D(segment_get_vertex_attrib(&prim, ivert, S2D_POSITION, &attr));
+ WRITE(attr.value, 2);
+ }
+ } else {
+ struct s3d_primitive prim;
+ S3D(scene_view_get_primitive(scn->s3d_view, (unsigned)iprim, &prim));
+ FOR_EACH(ivert, 0, 3) {
+ struct s3d_attrib attr;
+ S3D(triangle_get_vertex_attrib(&prim, ivert, S3D_POSITION, &attr));
+ WRITE(attr.value, 3);
+ }
+ }
+
+ interf = scene_get_interface(scn, (unsigned)iprim);
+ WRITE(&interf->medium_front->type, 1);
+ WRITE(&interf->medium_front->id, 1);
+ WRITE(&interf->medium_back->type, 1);
+ WRITE(&interf->medium_back->id, 1);
+ }
+ #undef WRITE
+
+ len = (size_t)ftell(stream);
+ rewind(stream);
+#ifdef COMPILER_GCC
+ data = mmap(NULL, len, PROT_READ, MAP_PRIVATE|MAP_POPULATE, fileno(stream), 0);
+ if(data == MAP_FAILED) {
+ res = RES_IO_ERR;
+ goto error;
+ }
+#else
+ data = MEM_ALLOC_ALIGNED(scn->dev->allocator, len, 16);
+ if(!data) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ if(fread(data, len, 1, stream) != 1) {
+ res = RES_IO_ERR;
+ goto error;
+ }
+#endif
+
+ res = hash_sha256(scn->dev->allocator, data, len, hash);
+ if(res != RES_OK) goto error;
+
+exit:
+#ifdef COMPILER_GCC
+ if(data) munmap(data, len);
+#else
+ if(data) MEM_RM(scn->dev->allocator, data);
+#endif
+ if(stream) fclose(stream);
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -24,6 +24,7 @@
#include "sdis_scene_c.h"
#include <star/ssp.h>
+#include <rsys/cstr.h>
#include <rsys/float22.h>
#include <rsys/float33.h>
#include <rsys/rsys.h>
@@ -957,6 +958,64 @@ error:
goto exit;
}
+static res_T
+XD(scene_find_closest_point)
+ (const struct sdis_scene* scn,
+ const double pos[3],
+ const double radius,
+ size_t* iprim,
+ double uv[2])
+{
+ struct sXd(hit) hit;
+ float query_pos[DIM];
+ float query_radius;
+ res_T res = RES_OK;
+
+ if(!scn || !pos || radius <= 0 || !iprim || !uv
+ || scene_is_2d(scn) != (DIM == 2)) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Avoid a null query radius due to casting in single-precision */
+ query_radius = MMAX((float)radius, FLT_MIN);
+
+ fX_set_dX(query_pos, pos);
+ res = sXd(scene_view_closest_point)
+ (scn->sXd(view), query_pos, query_radius, NULL, &hit);
+ if(res != RES_OK) {
+#if DIM == 2
+ log_err(scn->dev,
+ "%s: error querying the closest position at {%g, %g} "
+ "for a radius of %g -- %s.\n",
+ FUNC_NAME, SPLIT2(query_pos), query_radius, res_to_cstr(res));
+#else
+ log_err(scn->dev,
+ "%s: error querying the closest position at {%g, %g, %g} "
+ "for a radius of %g -- %s.\n",
+ FUNC_NAME, SPLIT3(query_pos), query_radius, res_to_cstr(res));
+#endif
+ goto error;
+ }
+
+ if(SXD_HIT_NONE(&hit)) {
+ *iprim = SDIS_PRIMITIVE_NONE;
+ } else {
+ *iprim = hit.prim.scene_prim_id;
+#if DIM == 2
+ uv[0] = hit.u;
+#else
+ uv[0] = hit.uv[0];
+ uv[1] = hit.uv[1];
+#endif
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -20,6 +20,7 @@
#include <star/s3d.h>
#include <rsys/dynamic_array_uint.h>
+#include <rsys/hash.h>
#include <rsys/hash_table.h>
#include <rsys/ref_count.h>
@@ -248,6 +249,11 @@ scene_get_medium_in_closed_boundaries
const double position[],
struct sdis_medium** medium);
+extern LOCAL_SYM res_T
+scene_compute_hash
+ (const struct sdis_scene* scn,
+ hash256_T hash);
+
static INLINE void
scene_get_enclosure_ids
(const struct sdis_scene* scn,
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -540,3 +540,16 @@ sdis_solve_medium_green_function
}
}
+res_T
+sdis_compute_power
+ (struct sdis_scene* scn,
+ const struct sdis_compute_power_args* args,
+ struct sdis_estimator** estimator)
+{
+ if(!scn) return RES_BAD_ARG;
+ if(scene_is_2d(scn)) {
+ return compute_power_2d(scn, args, estimator);
+ } else {
+ return compute_power_3d(scn, args, estimator);
+ }
+}
diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h
@@ -214,7 +214,7 @@ XD(solve_boundary)
greens = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*greens));
if(!greens) { res = RES_MEM_ERR; goto error; }
FOR_EACH(i, 0, scn->dev->nthreads) {
- res = green_function_create(scn->dev, &greens[i]);
+ res = green_function_create(scn, &greens[i]);
if(res != RES_OK) goto error;
}
}
diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h
@@ -284,7 +284,7 @@ XD(solve_medium)
greens = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*greens));
if(!greens) { res = RES_MEM_ERR; goto error; }
FOR_EACH(i, 0, scn->dev->nthreads) {
- res = green_function_create(scn->dev, &greens[i]);
+ res = green_function_create(scn, &greens[i]);
if(res != RES_OK) goto error;
}
}
@@ -477,4 +477,192 @@ error:
goto exit;
}
+static res_T
+XD(compute_power)
+ (struct sdis_scene* scn,
+ const struct sdis_compute_power_args* args,
+ struct sdis_estimator** out_estimator)
+{
+ struct darray_enclosure_cumul cumul;
+ struct sdis_estimator* estimator = NULL;
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ struct accum* acc_mpows = NULL;
+ struct accum* acc_times = NULL;
+ double spread = 0;
+ size_t i = 0;
+ size_t nrealisations = 0;
+ int64_t irealisation = 0;
+ int cumul_is_init = 0;
+ int progress = 0;
+ ATOMIC nsolved_realisations = 0;
+ ATOMIC res = RES_OK;
+
+ if(!scn
+ || !args
+ || !out_estimator
+ || !args->medium
+ || !args->nrealisations
+ || args->nrealisations > INT64_MAX
+ || args->fp_to_meter <= 0
+ || args->time_range[0] < 0
+ || args->time_range[0] > args->time_range[1]
+ || ( args->time_range[1] > DBL_MAX
+ && args->time_range[0] != args->time_range[1])) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(scene_is_2d(scn) != (SDIS_XD_DIMENSION==2)) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(sdis_medium_get_type(args->medium) != SDIS_SOLID) {
+ log_err(scn->dev, "Could not compute mean power on a non solid medium.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Create the RNG proxy */
+ if(args->rng_state) {
+ res = ssp_rng_proxy_create_from_rng(scn->dev->allocator, args->rng_state,
+ scn->dev->nthreads, &rng_proxy);
+ if(res != RES_OK) goto error;
+ } else {
+ res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
+ scn->dev->nthreads, &rng_proxy);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Create the per thread RNG */
+ rngs = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*rngs));
+ if(!rngs) { res = RES_MEM_ERR; goto error; }
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Create the per thread accumulators */
+ acc_mpows = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_mpows));
+ if(!acc_mpows) { res = RES_MEM_ERR; goto error; }
+ acc_times = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_times));
+ if(!acc_times) { res = RES_MEM_ERR; goto error; }
+
+ /* Compute the cumulative of the spreads of the enclosures surrounding the
+ * submitted medium */
+ darray_enclosure_cumul_init(scn->dev->allocator, &cumul);
+ cumul_is_init = 1;
+ res = compute_medium_enclosure_cumulative(scn, args->medium, &cumul);
+ if(res != RES_OK) goto error;
+
+ /* Fetch the overall medium spread from the unormalized enclosure cumulative */
+ spread = darray_enclosure_cumul_cdata_get(&cumul)
+ [darray_enclosure_cumul_size_get(&cumul)-1].cumul;
+
+ /* Create the estimator */
+ res = estimator_create(scn->dev, SDIS_ESTIMATOR_POWER, &estimator);
+ if(res != RES_OK) goto error;
+
+ nrealisations = args->nrealisations;
+ omp_set_num_threads((int)scn->dev->nthreads);
+ #pragma omp parallel for schedule(static)
+ for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ struct time t0, t1;
+ struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
+ const int ithread = omp_get_thread_num();
+ struct ssp_rng* rng = rngs[ithread];
+ struct accum* acc_mpow = &acc_mpows[ithread];
+ struct accum* acc_time = &acc_times[ithread];
+ const struct enclosure* enc = NULL;
+ double power = 0;
+ double usec = 0;
+ size_t n = 0;
+ int pcent = 0;
+ res_T res_local = RES_OK;
+
+ if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+
+ /* Begin time registration */
+ time_current(&t0);
+
+ /* Sample the time */
+ vtx.time = sample_time(rng, args->time_range);
+
+ /* Uniformly Sample an enclosure that surround the submitted medium and
+ * uniformly sample a position into it */
+ enc = sample_medium_enclosure(&cumul, rng);
+ res_local = XD(sample_enclosure_position)(enc, rng, vtx.P);
+ if(res_local != RES_OK) {
+ log_err(scn->dev, "%s: could not sample a medium position.\n", FUNC_NAME);
+ ATOMIC_SET(&res, res_local);
+ goto error_it;
+ }
+
+ /* Fetch the volumic power */
+ power = solid_get_volumic_power(args->medium, &vtx);
+
+ /* Stop time registration */
+ time_sub(&t0, time_current(&t1), &t0);
+ usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+
+ /* Update accumulators */
+ acc_mpow->sum += power; acc_mpow->sum2 += power*power; ++acc_mpow->count;
+ acc_time->sum += usec; acc_time->sum2 += usec*usec; ++acc_time->count;
+
+ /* Update progress */
+ n = (size_t)ATOMIC_INCR(&nsolved_realisations);
+ pcent = (int)((double)n * 100.0 / (double)nrealisations + 0.5/*round*/);
+ #pragma omp critical
+ if(pcent > progress) {
+ progress = pcent;
+ log_info(scn->dev, "Computing mean power: %3d%%\r", progress);
+ }
+ exit_it:
+ continue;
+ error_it:
+ goto exit_it;
+ }
+ if(res != RES_OK) goto error;
+
+ /* Add a new line after the progress status */
+ log_info(scn->dev, "Computing mean power: %3d%%\n", progress);
+
+ /* Setup the estimated mean power */
+ {
+ struct accum acc_mpow;
+ struct accum acc_time;
+ sum_accums(acc_mpows, scn->dev->nthreads, &acc_mpow);
+ sum_accums(acc_times, scn->dev->nthreads, &acc_time);
+ ASSERT(acc_mpow.count == acc_time.count);
+
+ estimator_setup_realisations_count(estimator, nrealisations, acc_mpow.count);
+ estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2);
+ estimator_setup_power
+ (estimator, acc_mpow.sum, acc_mpow.sum2, spread, args->time_range);
+ res = estimator_save_rng_state(estimator, rng_proxy);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ if(rngs) {
+ FOR_EACH(i, 0, scn->dev->nthreads) {if(rngs[i]) SSP(rng_ref_put(rngs[i]));}
+ MEM_RM(scn->dev->allocator, rngs);
+ }
+ if(acc_mpows) MEM_RM(scn->dev->allocator, acc_mpows);
+ if(acc_times) MEM_RM(scn->dev->allocator, acc_times);
+ if(cumul_is_init) darray_enclosure_cumul_release(&cumul);
+ if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ if(out_estimator) *out_estimator = estimator;
+ return (res_T)res;
+error:
+ if(estimator) {
+ SDIS(estimator_ref_put(estimator));
+ estimator = NULL;
+ }
+ goto exit;
+}
+
#include "sdis_Xd_end.h"
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -115,7 +115,7 @@ XD(solve_probe)
greens = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*greens));
if(!greens) { res = RES_MEM_ERR; goto error; }
FOR_EACH(i, 0, scn->dev->nthreads) {
- res = green_function_create(scn->dev, &greens[i]);
+ res = green_function_create(scn, &greens[i]);
if(res != RES_OK) goto error;
}
}
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -150,7 +150,7 @@ XD(solve_probe_boundary)
greens = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*greens));
if(!greens) { res = RES_MEM_ERR; goto error; }
FOR_EACH(i, 0, scn->dev->nthreads) {
- res = green_function_create(scn->dev, &greens[i]);
+ res = green_function_create(scn, &greens[i]);
if(res != RES_OK) goto error;
}
}
diff --git a/src/test_sdis_compute_mean_power.c b/src/test_sdis_compute_mean_power.c
@@ -0,0 +1,343 @@
+/* Copyright (C) 2016-2020 |Meso|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/stretchy_array.h>
+#include <star/s3dut.h>
+
+#define UNKOWN_TEMPERATURE -1
+#define N 100000ul /* #realisations */
+#define POWER0 10
+#define POWER1 5
+
+static INLINE void
+check_intersection
+ (const double val0,
+ const double eps0,
+ const double val1,
+ const double eps1)
+{
+ double interval0[2], interval1[2];
+ double intersection[2];
+ interval0[0] = val0 - eps0;
+ interval0[1] = val0 + eps0;
+ interval1[0] = val1 - eps1;
+ interval1[1] = val1 + eps1;
+ intersection[0] = MMAX(interval0[0], interval1[0]);
+ intersection[1] = MMIN(interval0[1], interval1[1]);
+ CHK(intersection[0] <= intersection[1]);
+}
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct context {
+ struct s3dut_mesh_data msh0;
+ struct s3dut_mesh_data msh1;
+ struct sdis_interface* interf0;
+ struct sdis_interface* interf1;
+};
+
+static void
+get_indices(const size_t itri, size_t ids[3], void* context)
+{
+ const struct context* ctx = context;
+ /* Note that we swap the indices to ensure that the triangle normals point
+ * inward the super shape */
+ if(itri < ctx->msh0.nprimitives) {
+ ids[0] = ctx->msh0.indices[itri*3+0];
+ ids[2] = ctx->msh0.indices[itri*3+1];
+ ids[1] = ctx->msh0.indices[itri*3+2];
+ } else {
+ const size_t itri2 = itri - ctx->msh0.nprimitives;
+ ids[0] = ctx->msh1.indices[itri2*3+0] + ctx->msh0.nvertices;
+ ids[2] = ctx->msh1.indices[itri2*3+1] + ctx->msh0.nvertices;
+ ids[1] = ctx->msh1.indices[itri2*3+2] + ctx->msh0.nvertices;
+ }
+}
+
+static void
+get_position(const size_t ivert, double pos[3], void* context)
+{
+ const struct context* ctx = context;
+ if(ivert < ctx->msh0.nvertices) {
+ pos[0] = ctx->msh0.positions[ivert*3+0] - 2.0;
+ pos[1] = ctx->msh0.positions[ivert*3+1];
+ pos[2] = ctx->msh0.positions[ivert*3+2];
+ } else {
+ const size_t ivert2 = ivert - ctx->msh0.nvertices;
+ pos[0] = ctx->msh1.positions[ivert2*3+0] + 2.0;
+ pos[1] = ctx->msh1.positions[ivert2*3+1];
+ pos[2] = ctx->msh1.positions[ivert2*3+2];
+ }
+}
+
+static void
+get_interface(const size_t itri, struct sdis_interface** bound, void* context)
+{
+ const struct context* ctx = context;
+ *bound = itri < ctx->msh0.nprimitives ? ctx->interf0 : ctx->interf1;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag != NULL); (void)data;
+ return 1.0;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx == NULL); (void)data;
+ return 300;
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 1.0;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 1.0;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 1.0;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 1.0 / 20.0;
+}
+
+static double
+solid_get_volumic_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return vtx->P[0] < 0 ? POWER0 : POWER1;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct context ctx;
+ struct s3dut_mesh* sphere = NULL;
+ struct s3dut_mesh* cylinder = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid0 = NULL;
+ struct sdis_medium* solid1 = NULL;
+ struct sdis_interface* interf0 = NULL;
+ struct sdis_interface* interf1 = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_mc mpow = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_compute_power_args args = SDIS_COMPUTE_POWER_ARGS_DEFAULT;
+ size_t nverts = 0;
+ size_t ntris = 0;
+ double ref = 0;
+ (void)argc, (void) argv;
+
+ OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+
+ /* Setup the fluid shader */
+ fluid_shader.temperature = fluid_get_temperature;
+
+ /* 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 = solid_get_delta;
+ solid_shader.volumic_power = solid_get_volumic_power;
+
+ /* Create the fluid */
+ OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
+
+ /* Create the solids */
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid0));
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid1));
+
+ /* Create the interface0 */
+ OK(sdis_interface_create(dev, solid0, fluid, &interf_shader, NULL, &interf0));
+ OK(sdis_interface_create(dev, solid1, fluid, &interf_shader, NULL, &interf1));
+ ctx.interf0 = interf0;
+ ctx.interf1 = interf1;
+
+ /* Create the geometry */
+ OK(s3dut_create_sphere(&allocator, 1, 512, 256, &sphere));
+ OK(s3dut_create_cylinder(&allocator, 1, 10, 512, 8, &cylinder));
+ OK(s3dut_mesh_get_data(sphere, &ctx.msh0));
+ OK(s3dut_mesh_get_data(cylinder, &ctx.msh1));
+
+ /* Create the scene */
+ ntris = ctx.msh0.nprimitives + ctx.msh1.nprimitives;
+ nverts = ctx.msh0.nvertices + ctx.msh1.nvertices;
+ OK(sdis_scene_create(dev, ntris, get_indices, get_interface, nverts,
+ get_position, &ctx, &scn));
+
+ /* Test sdis_compute_power function */
+ args.nrealisations = N;
+ args.medium = solid0;
+ args.time_range[0] = INF;
+ args.time_range[1] = INF;
+ BA(sdis_compute_power(NULL, &args, &estimator));
+ BA(sdis_compute_power(scn, NULL, &estimator));
+ BA(sdis_compute_power(scn, &args, NULL));
+ args.nrealisations = 0;
+ BA(sdis_compute_power(scn, &args, &estimator));
+ args.nrealisations = N;
+ args.medium = NULL;
+ BA(sdis_compute_power(scn, &args, &estimator));
+ args.medium = solid0;
+ args.fp_to_meter = 0;
+ BA(sdis_compute_power(scn, &args, &estimator));
+ args.fp_to_meter = 1;
+ args.time_range[0] = args.time_range[1] = -1;
+ BA(sdis_compute_power(scn, &args, &estimator));
+ args.time_range[0] = 1;
+ BA(sdis_compute_power(scn, &args, &estimator));
+ args.time_range[1] = 0;
+ BA(sdis_compute_power(scn, &args, &estimator));
+ args.time_range[0] = args.time_range[1] = INF;
+ OK(sdis_compute_power(scn, &args, &estimator));
+
+ BA(sdis_estimator_get_power(NULL, &mpow));
+ BA(sdis_estimator_get_power(estimator, NULL));
+ OK(sdis_estimator_get_power(estimator, &mpow));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
+ /* Check results for solid 0 */
+ ref = 4.0/3.0 * PI * POWER0;
+ printf("Mean power of the solid0 = %g ~ %g +/- %g\n",
+ ref, mpow.E, mpow.SE);
+ check_intersection(ref, 1.e-2, mpow.E, 3*mpow.SE);
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Check results for solid 1 */
+ args.medium = solid1;
+ OK(sdis_compute_power(scn, &args, &estimator));
+ OK(sdis_estimator_get_power(estimator, &mpow));
+ ref = PI * 10 * POWER1;
+ printf("Mean power of the solid1 = %g ~ %g +/- %g\n",
+ ref, mpow.E, mpow.SE);
+ check_intersection(ref, 1.e-2, mpow.E, 3*mpow.SE);
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Check for a not null time range */
+ args.time_range[0] = 0;
+ args.time_range[1] = 10;
+ OK(sdis_compute_power(scn, &args, &estimator));
+ OK(sdis_estimator_get_power(estimator, &mpow));
+ ref = PI * 10 * POWER1 / 10;
+ printf("Mean power of the solid1 in [0, 10] s = %g ~ %g +/- %g\n",
+ ref, mpow.E, mpow.SE);
+ check_intersection(ref, 1.e-2, mpow.E, 3*mpow.SE);
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Reset the scene with only one solid medium */
+ OK(sdis_scene_ref_put(scn));
+ ctx.interf0 = interf0;
+ ctx.interf1 = interf0;
+ OK(sdis_scene_create(dev, ntris, get_indices, get_interface, nverts,
+ get_position, &ctx, &scn));
+
+ /* Check invalid medium */
+ args.time_range[0] = args.time_range[1] = 1;
+ args.medium = solid1;
+ BA(sdis_compute_power(scn, &args, &estimator));
+
+ /* Check non constant volumic power */
+ args.medium = solid0;
+ OK(sdis_compute_power(scn, &args, &estimator));
+ OK(sdis_estimator_get_power(estimator, &mpow));
+ ref = 4.0/3.0*PI*POWER0 + PI*10*POWER1;
+ printf("Mean power of the sphere+cylinder = %g ~ %g +/- %g\n",
+ ref, mpow.E, mpow.SE);
+ check_intersection(ref, 1.5e-1, mpow.E, 3*mpow.SE);
+ OK(sdis_estimator_ref_put(estimator));
+
+#if 0
+ {
+ double* vertices = NULL;
+ size_t* indices = NULL;
+ size_t i;
+ CHK(vertices = MEM_CALLOC(&allocator, nverts*3, sizeof(*vertices)));
+ CHK(indices = MEM_CALLOC(&allocator, ntris*3, sizeof(*indices)));
+ FOR_EACH(i, 0, ntris) get_indices(i, indices + i*3, &ctx);
+ FOR_EACH(i, 0, nverts) get_position(i, vertices + i*3, &ctx);
+ dump_mesh(stdout, vertices, nverts, indices, ntris);
+ MEM_RM(&allocator, vertices);
+ MEM_RM(&allocator, indices);
+ }
+#endif
+
+ /* Clean up memory */
+ OK(sdis_device_ref_put(dev));
+ OK(sdis_medium_ref_put(fluid));
+ OK(sdis_medium_ref_put(solid0));
+ OK(sdis_medium_ref_put(solid1));
+ OK(sdis_interface_ref_put(interf0));
+ OK(sdis_interface_ref_put(interf1));
+ OK(sdis_scene_ref_put(scn));
+ OK(s3dut_mesh_ref_put(sphere));
+ OK(s3dut_mesh_ref_put(cylinder));
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -278,6 +278,7 @@ main(int argc, char** argv)
const double T1 = 310; /* Fixed temperature on the right side of the system */
const double thickness = 2.0; /* Thickness of the solid along X */
double Ts0, Ts1, hr, tmp;
+ struct interfac* p_intface;
(void)argc, (void)argv;
OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
@@ -376,11 +377,10 @@ main(int argc, char** argv)
get_position, &geom, &scn));
hr = 4.0 * BOLTZMANN_CONSTANT * Tref*Tref*Tref * emissivity;
- tmp = lambda/(2*lambda + thickness*hr) * (T1 - T0);
- Ts0 = T0 + tmp;
- Ts1 = T1 - tmp;
/* Run the simulations */
+ p_intface
+ = (struct interfac*)sdis_data_get(sdis_interface_get_data(interfaces[4]));
OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng));
FOR_EACH(isimul, 0, nsimuls) {
struct sdis_mc T = SDIS_MC_NULL;
@@ -393,6 +393,10 @@ main(int argc, char** argv)
size_t nreals = 0;
size_t nfails = 0;
const size_t N = 10000;
+ double T1b;
+
+ /* Reset temperature */
+ p_intface->temperature = T1;
solve_args.nrealisations = N;
solve_args.time_range[0] = INF;
@@ -400,6 +404,7 @@ main(int argc, char** argv)
solve_args.position[0] = ssp_rng_uniform_double(rng, -0.9, 0.9);
solve_args.position[1] = ssp_rng_uniform_double(rng, -0.9, 0.9);
solve_args.position[2] = ssp_rng_uniform_double(rng, -0.9, 0.9);
+ u = (solve_args.position[0] + 1) / thickness;
solve_args.reference_temperature = Tref;
OK(sdis_solve_probe(scn, &solve_args, &estimator));
@@ -408,10 +413,12 @@ main(int argc, char** argv)
OK(sdis_estimator_get_temperature(estimator, &T));
OK(sdis_estimator_get_realisation_time(estimator, &time));
- u = (solve_args.position[0] + 1) / thickness;
+ tmp = lambda / (2 * lambda + thickness * hr) * (T1 - T0);
+ Ts0 = T0 + tmp;
+ Ts1 = T1 - tmp;
ref = u * Ts1 + (1-u) * Ts0;
- printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
- SPLIT3(solve_args.position), ref, T.E, T.SE);
+ printf("Temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n",
+ SPLIT3(solve_args.position), p_intface->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);
@@ -424,6 +431,37 @@ main(int argc, char** argv)
OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, scn, solve_args.time_range);
+
+ OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
+ printf("\n");
+
+ /* Check green used at a different temperature */
+ p_intface->temperature = T1b = T1 + ((double)isimul + 1) * 10;
+
+ OK(sdis_solve_probe(scn, &solve_args, &estimator));
+ 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));
+
+ tmp = lambda / (2 * lambda + thickness * hr) * (T1b - T0);
+ Ts0 = T0 + tmp;
+ Ts1 = T1b - tmp;
+ ref = u * Ts1 + (1 - u) * Ts0;
+ printf("Temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n",
+ SPLIT3(solve_args.position), p_intface->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);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails < N / 1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE) == 1);
+
+ OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
+ check_green_function(green);
+ check_estimator_eq(estimator, estimator2);
OK(sdis_estimator_ref_put(estimator));
OK(sdis_estimator_ref_put(estimator2));
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -415,7 +415,7 @@ main(int argc, char** argv)
u = (solve_args.position[0] + 1) / thickness;
ref = u * Ts1 + (1-u) * Ts0;
- printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
+ printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
SPLIT2(solve_args.position), 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);
@@ -429,6 +429,7 @@ main(int argc, char** argv)
OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, scn, solve_args.time_range);
OK(sdis_estimator_ref_put(estimator));
OK(sdis_estimator_ref_put(estimator2));
diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c
@@ -300,6 +300,7 @@ main(int argc, char** argv)
OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, box_scn, solve_args.time_range);
OK(sdis_estimator_ref_put(estimator2));
OK(sdis_green_function_ref_put(green));
}
@@ -341,6 +342,7 @@ main(int argc, char** argv)
OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, square_scn, solve_args.time_range);
OK(sdis_estimator_ref_put(estimator2));
OK(sdis_green_function_ref_put(green));
}
diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c
@@ -316,6 +316,7 @@ main(int argc, char** argv)
OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, box_scn, solve_args.time_range);
OK(sdis_estimator_ref_put(estimator2));
OK(sdis_green_function_ref_put(green));
}
@@ -356,6 +357,7 @@ main(int argc, char** argv)
OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, square_scn, solve_args.time_range);
OK(sdis_estimator_ref_put(estimator2));
OK(sdis_green_function_ref_put(green));
}
diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c
@@ -127,7 +127,10 @@ interface_get_flux
* Helper functions
******************************************************************************/
static void
-solve(struct sdis_scene* scn, const double pos[])
+solve
+ (struct sdis_scene* scn,
+ const double pos[],
+ struct interf* interf)
{
char dump[128];
struct time t0, t1, t2;
@@ -142,9 +145,12 @@ solve(struct sdis_scene* scn, const double pos[])
double ref;
const double time_range[2] = {INF, INF};
enum sdis_scene_dimension dim;
- ASSERT(scn && pos);
+ ASSERT(scn && pos && interf);
- ref = T0 + (1 - pos[0]) * PHI/LAMBDA;
+ /* Restore phi value */
+ interf->phi = PHI;
+
+ ref = T0 + (1 - pos[0]) * interf->phi/LAMBDA;
OK(sdis_scene_get_dimension(scn, &dim));
@@ -167,18 +173,18 @@ solve(struct sdis_scene* scn, const double pos[])
switch(dim) {
case SDIS_SCENE_2D:
- printf("Temperature at (%g %g) = %g ~ %g +/- %g\n",
- SPLIT2(pos), ref, T.E, T.SE);
+ printf("Temperature at (%g %g) with phi=%g = %g ~ %g +/- %g\n",
+ SPLIT2(pos), interf->phi, ref, T.E, T.SE);
break;
case SDIS_SCENE_3D:
- printf("Temperature at (%g %g %g) = %g ~ %g +/- %g\n",
- SPLIT3(pos), ref, T.E, T.SE);
+ printf("Temperature at (%g %g %g) with phi=%g = %g ~ %g +/- %g\n",
+ SPLIT3(pos), interf->phi, ref, T.E, T.SE);
break;
default: FATAL("Unreachable code.\n"); break;
}
- printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- printf("Elapsed time = %s\n\n", dump);
+ 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);
@@ -196,12 +202,12 @@ solve(struct sdis_scene* scn, const double pos[])
switch(dim) {
case SDIS_SCENE_2D:
- printf("Green temperature at (%g %g) = %g ~ %g +/- %g\n",
- SPLIT2(pos), ref, T.E, T.SE);
+ printf("Green temperature at (%g %g) with phi=%g = %g ~ %g +/- %g\n",
+ SPLIT2(pos), interf->phi, ref, T.E, T.SE);
break;
case SDIS_SCENE_3D:
- printf("Green temperature at (%g %g %g) = %g ~ %g +/- %g\n",
- SPLIT3(pos), ref, T.E, T.SE);
+ printf("Green temperature at (%g %g %g) with phi=%g = %g ~ %g +/- %g\n",
+ SPLIT3(pos), interf->phi, ref, T.E, T.SE);
break;
default: FATAL("Unreachable code.\n"); break;
}
@@ -215,6 +221,54 @@ solve(struct sdis_scene* scn, const double pos[])
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, scn, time_range);
+
+ OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
+ printf("\n");
+
+ /* Check green used at a different phi */
+ interf->phi = 3 * PHI;
+
+ 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));
+
+ ref = T0 + (1 - pos[0]) * interf->phi/LAMBDA;
+
+ switch (dim) {
+ case SDIS_SCENE_2D:
+ printf("Temperature at (%g %g) with phi=%g = %g ~ %g +/- %g\n",
+ SPLIT2(pos), interf->phi, ref, T.E, T.SE);
+ break;
+ case SDIS_SCENE_3D:
+ printf("Temperature at (%g %g %g) with phi=%g = %g ~ %g +/- %g\n",
+ SPLIT3(pos), interf->phi, ref, 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, T.SE * 3));
+
+ time_current(&t0);
+ OK(sdis_green_function_solve(green, time_range, &estimator2));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+ printf("Green solve time = %s\n", dump);
+
+ check_green_function(green);
+ check_estimator_eq(estimator, estimator2);
OK(sdis_estimator_ref_put(estimator));
OK(sdis_estimator_ref_put(estimator2));
@@ -328,9 +382,9 @@ main(int argc, char** argv)
/* Solve */
d3_splat(pos, 0.25);
printf(">> Box scene\n");
- solve(box_scn, pos);
+ solve(box_scn, pos, interf_props);
printf(">> Square Scene\n");
- solve(square_scn, pos);
+ solve(square_scn, pos, interf_props);
OK(sdis_scene_ref_put(box_scn));
OK(sdis_scene_ref_put(square_scn));
diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c
@@ -89,17 +89,19 @@ 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)
{
- struct sdis_scene* scn = NULL;
+ size_t duplicated_indices[] = { 0, 1, 2, 0, 2, 1 };
+ size_t degenerated_indices[] = { 0, 1, 1 };
+ double duplicated_vertices[] = { 0, 0, 0, 1, 1, 1, 0, 0, 0 };
double lower[3], upper[3];
- double uv0[2], uv1[2], pos[3], pos1[3];
+ double uv0[2], uv1[2], uv2[2], pos[3], pos1[3];
struct context ctx;
struct senc2d_scene* scn2d;
struct senc3d_scene* scn3d;
+ struct sdis_scene* scn = NULL;
size_t ntris, npos;
+ size_t iprim;
size_t i;
- size_t duplicated_indices[] = { 0, 1, 2, 0, 2, 1 };
- size_t degenerated_indices[] = { 0, 1, 1 };
- double duplicated_vertices[] = { 0, 0, 0, 1, 1, 1, 0, 0, 0 };
+ double dst = 0;
size_t dup_vrtx_indices[] = { 0, 1, 2 };
enum sdis_scene_dimension dim;
@@ -172,7 +174,32 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
BA(sdis_scene_boundary_project_position(scn, 6, pos, NULL));
OK(sdis_scene_boundary_project_position(scn, 6, pos, uv1));
+ BA(sdis_scene_find_closest_point(NULL, pos, INF, &iprim, uv2));
+ BA(sdis_scene_find_closest_point(scn, NULL, INF, &iprim, uv2));
+ BA(sdis_scene_find_closest_point(scn, pos, 0, &iprim, uv2));
+ BA(sdis_scene_find_closest_point(scn, pos, INF, NULL, uv2));
+ BA(sdis_scene_find_closest_point(scn, pos, INF, &iprim, NULL));
+ OK(sdis_scene_find_closest_point(scn, pos, INF, &iprim, uv2));
+
+ CHK(iprim == 6);
CHK(d2_eq_eps(uv0, uv1, 1.e-6));
+ CHK(d2_eq_eps(uv1, uv2, 1.e-6));
+
+ pos[0] = 0.5;
+ pos[1] = 0.1;
+ pos[2] = 0.25;
+ OK(sdis_scene_find_closest_point(scn, pos, INF, &iprim, uv2));
+ CHK(iprim == 10);
+
+ OK(sdis_scene_boundary_project_position(scn, 10, pos, uv0));
+ CHK(d2_eq_eps(uv0, uv2, 1.e-6));
+
+ OK(sdis_scene_get_boundary_position(scn, iprim, uv2, pos1));
+ dst = d3_len(d3_sub(pos1, pos, pos1));
+ CHK(eq_eps(dst, 0.1, 1.e-6));
+
+ OK(sdis_scene_find_closest_point(scn, pos, 0.09, &iprim, uv2));
+ CHK(iprim == SDIS_PRIMITIVE_NONE);
FOR_EACH(i, 0, 64) {
uv0[0] = rand_canonic();
@@ -180,7 +207,10 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
OK(sdis_scene_get_boundary_position(scn, 4, uv0, pos));
OK(sdis_scene_boundary_project_position(scn, 4, pos, uv1));
+ OK(sdis_scene_find_closest_point(scn, pos, INF, &iprim, uv2));
CHK(d2_eq_eps(uv0, uv1, 1.e-6));
+ CHK(d2_eq_eps(uv1, uv2, 1.e-6));
+ CHK(iprim == 4);
}
pos[0] = 10;
@@ -208,17 +238,19 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
static void
test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
{
+ size_t duplicated_indices[] = { 0, 1, 1, 0 };
+ size_t degenerated_indices[] = { 0, 0 };
+ double duplicated_vertices[] = { 0, 0, 0, 0 };
struct sdis_scene* scn = NULL;
double lower[2], upper[2];
- double u0, u1, pos[2];
+ double u0, u1, u2, pos[2], pos1[2];
+ double dst;
struct context ctx;
struct senc2d_scene* scn2d;
struct senc3d_scene* scn3d;
size_t nsegs, npos;
size_t i;
- size_t duplicated_indices[] = { 0, 1, 1, 0 };
- size_t degenerated_indices[] = { 0, 0 };
- double duplicated_vertices[] = { 0, 0, 0, 0 };
+ size_t iprim;
size_t dup_vrtx_indices[] = { 0, 1 };
enum sdis_scene_dimension dim;
@@ -288,14 +320,41 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
BA(sdis_scene_boundary_project_position(scn, 1, pos, NULL));
OK(sdis_scene_boundary_project_position(scn, 1, pos, &u1));
+ BA(sdis_scene_find_closest_point(NULL, pos, INF, &iprim, &u2));
+ BA(sdis_scene_find_closest_point(scn, NULL, INF, &iprim, &u2));
+ BA(sdis_scene_find_closest_point(scn, pos, 0, &iprim, &u2));
+ BA(sdis_scene_find_closest_point(scn, pos, INF, NULL, &u2));
+ BA(sdis_scene_find_closest_point(scn, pos, INF, &iprim, NULL));
+ OK(sdis_scene_find_closest_point(scn, pos, INF, &iprim, &u2));
+
CHK(eq_eps(u0, u1, 1.e-6));
+ CHK(eq_eps(u1, u2, 1.e-6));
+ CHK(iprim == 1);
+
+ pos[0] = 0.5;
+ pos[1] = 0.1;
+ OK(sdis_scene_find_closest_point(scn, pos, INF, &iprim, &u2));
+ CHK(iprim == 0);
+
+ OK(sdis_scene_boundary_project_position(scn, 0, pos, &u0));
+ CHK(eq_eps(u0, u2, 1.e-6));
+
+ OK(sdis_scene_get_boundary_position(scn, iprim, &u2, pos1));
+ dst = d2_len(d2_sub(pos1, pos, pos1));
+ CHK(eq_eps(dst, 0.1, 1.e-6));
+
+ OK(sdis_scene_find_closest_point(scn, pos, 0.09, &iprim, &u2));
+ CHK(iprim == SDIS_PRIMITIVE_NONE);
FOR_EACH(i, 0, 64) {
u0 = rand_canonic();
OK(sdis_scene_get_boundary_position(scn, 2, &u0, pos));
OK(sdis_scene_boundary_project_position(scn, 2, pos, &u1));
+ OK(sdis_scene_find_closest_point(scn, pos, INF, &iprim, &u2));
CHK(eq_eps(u0, u1, 1.e-6));
+ CHK(eq_eps(u1, u2, 1.e-6));
+ CHK(iprim == 2);
}
d2(pos, 5, 0.5);
diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c
@@ -348,6 +348,7 @@ main(int argc, char** argv)
check_green_function(green);
OK(sdis_green_function_solve(green, probe_args.time_range, &estimator2));
check_estimator(estimator2, N, ref);
+ check_green_serialization(green, box_scn, probe_args.time_range);
OK(sdis_green_function_ref_put(green));
OK(sdis_estimator_ref_put(estimator));
@@ -378,6 +379,7 @@ main(int argc, char** argv)
check_green_function(green);
OK(sdis_green_function_solve(green, probe_args.time_range, &estimator2));
check_estimator(estimator2, N, ref);
+ check_green_serialization(green, square_scn, probe_args.time_range);
OK(sdis_estimator_ref_put(estimator));
OK(sdis_estimator_ref_put(estimator2));
@@ -465,11 +467,11 @@ main(int argc, char** argv)
BA(GREEN(box_scn, &bound_args, &green));
sides[0] = SDIS_FRONT;
-
OK(GREEN(box_scn, &bound_args, &green));
check_green_function(green);
OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2));
check_estimator(estimator2, N, ref);
+ check_green_serialization(green, box_scn, bound_args.time_range);
OK(sdis_green_function_ref_put(green));
OK(sdis_estimator_ref_put(estimator));
@@ -506,6 +508,7 @@ main(int argc, char** argv)
check_green_function(green);
OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2));
check_estimator(estimator2, N, ref);
+ check_green_serialization(green, square_scn, bound_args.time_range);
OK(sdis_green_function_ref_put(green));
OK(sdis_estimator_ref_put(estimator));
@@ -538,6 +541,7 @@ main(int argc, char** argv)
check_green_function(green);
OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2));
check_estimator(estimator2, N, ref);
+ check_green_serialization(green, box_scn, bound_args.time_range);
OK(sdis_green_function_ref_put(green));
OK(sdis_estimator_ref_put(estimator));
@@ -555,6 +559,7 @@ main(int argc, char** argv)
check_green_function(green);
OK(sdis_green_function_solve(green, bound_args.time_range, &estimator2));
check_estimator(estimator2, N, ref);
+ check_green_serialization(green, square_scn, bound_args.time_range);
OK(sdis_green_function_ref_put(green));
OK(sdis_estimator_ref_put(estimator));
diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c
@@ -459,6 +459,7 @@ main(int argc, char** argv)
OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, scn, solve_args.time_range);
OK(sdis_green_function_ref_put(green));
diff --git a/src/test_sdis_solve_medium_2d.c b/src/test_sdis_solve_medium_2d.c
@@ -402,6 +402,7 @@ main(int argc, char** argv)
OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, scn, solve_args.time_range);
OK(sdis_green_function_ref_put(green));
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -259,6 +259,7 @@ main(int argc, char** argv)
struct sdis_data* data = NULL;
struct sdis_estimator* estimator = NULL;
struct sdis_estimator* estimator2 = NULL;
+ struct sdis_estimator* estimator3 = NULL;
struct sdis_green_function* green = NULL;
const struct sdis_heat_path* path = NULL;
struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
@@ -272,6 +273,7 @@ main(int argc, char** argv)
struct interf* interface_param;
struct ssp_rng* rng_state = NULL;
enum sdis_estimator_type type;
+ FILE* stream = NULL;
double ref;
const size_t N = 1000;
const size_t N_dump = 10;
@@ -403,8 +405,8 @@ main(int argc, char** argv)
OK(sdis_estimator_get_rng_state(estimator, &rng_state));
ref = 300;
- printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
- SPLIT3(solve_args.position), ref, T.E, T.SE);
+ 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);
printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
@@ -444,14 +446,60 @@ main(int argc, char** argv)
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
+ printf("\n");
+
+ /* Check green used at a different temperature */
+ fluid_param->temperature = 500;
+ OK(sdis_solve_probe(scn, &solve_args, &estimator));
+ 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));
+
+ 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);
+ printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails < N / 1000);
+ CHK(eq_eps(T.E, ref, T.SE));
+
+ OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
+ check_green_function(green);
+ check_estimator_eq(estimator, estimator2);
+
+ CHK(stream = tmpfile());
+ BA(sdis_green_function_write(NULL, stream));
+ BA(sdis_green_function_write(green, NULL));
+ OK(sdis_green_function_write(green, stream));
+
BA(sdis_green_function_ref_get(NULL));
OK(sdis_green_function_ref_get(green));
BA(sdis_green_function_ref_put(NULL));
OK(sdis_green_function_ref_put(green));
OK(sdis_green_function_ref_put(green));
+ rewind(stream);
+ BA(sdis_green_function_create_from_stream(NULL, stream, &green));
+ BA(sdis_green_function_create_from_stream(scn, NULL, &green));
+ BA(sdis_green_function_create_from_stream(scn, stream, NULL));
+ OK(sdis_green_function_create_from_stream(scn, stream, &green));
+ CHK(!fclose(stream));
+
+ OK(sdis_green_function_solve(green, solve_args.time_range, &estimator3));
+
+ check_green_function(green);
+ check_estimator_eq_strict(estimator2, estimator3);
+
+ OK(sdis_green_function_ref_put(green));
+
OK(sdis_estimator_ref_put(estimator));
OK(sdis_estimator_ref_put(estimator2));
+ OK(sdis_estimator_ref_put(estimator3));
OK(sdis_solve_probe(scn, &solve_args, &estimator));
BA(sdis_estimator_get_paths_count(NULL, &n));
diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c
@@ -218,9 +218,9 @@ main(int argc, char** argv)
/* Setup the per primitive scene interfaces */
CHK(sizeof(interfaces)/sizeof(struct sdis_interface*) == box_ntriangles);
- interfaces[0] = interfaces[1] = T300; /* Front face */
+ interfaces[0] = interfaces[1] = T300; /* Back face */
interfaces[2] = interfaces[3] = Tnone; /* Left face */
- interfaces[4] = interfaces[5] = T350; /* Back face */
+ interfaces[4] = interfaces[5] = T350; /* Front face */
interfaces[6] = interfaces[7] = Tnone; /* Right face */
interfaces[8] = interfaces[9] = Tnone; /* Top face */
interfaces[10] = interfaces[11] = Tnone; /* Bottom face */
@@ -268,6 +268,7 @@ main(int argc, char** argv)
OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, scn, solve_args.time_range);
/* Release data */
OK(sdis_estimator_ref_put(estimator));
diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c
@@ -267,6 +267,7 @@ main(int argc, char** argv)
OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, scn, solve_args.time_range);
/* Release data */
OK(sdis_estimator_ref_put(estimator));
@@ -280,3 +281,4 @@ main(int argc, char** argv)
CHK(mem_allocated_size() == 0);
return 0;
}
+
diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c
@@ -324,6 +324,7 @@ main(int argc, char** argv)
OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, scn, solve_args.time_range);
/* Release data */
OK(sdis_estimator_ref_put(estimator));
diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c
@@ -316,6 +316,7 @@ main(int argc, char** argv)
OK(sdis_green_function_solve(green, solve_args.time_range, &estimator2));
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, scn, solve_args.time_range);
/* Release data */
OK(sdis_estimator_ref_put(estimator));
diff --git a/src/test_sdis_transcient.c b/src/test_sdis_transcient.c
@@ -0,0 +1,676 @@
+/* Copyright (C) 2016-2019 |Meso|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 <string.h>
+
+#define UNKNOWN_TEMPERATURE -1
+
+/*
+ * The scene is composed of a solid cuboid whose temperature is fixed on its 6
+ * faces. The test consist in checking that the estimated temperature at a
+ * given temperature and time is compatible with the reference temperature
+ * computed by analytically evaluating the green function.
+ *
+ * The test is performed on 3 scenes that actually represent the same system.
+ * The first scene is simply the cuboid, as it. The second scene is the same
+ * cuboid but this time formed by 2 sub cuboid with strictly the same physical
+ * properties. Finally, the last scene is the same cube with successive
+ * cubes into it used only to add fictive boundaries.
+ */
+
+static const double vertices[12/*#vertices*/*3/*#coords per vertex*/] = {
+ 0.0, 0.0, 0.0,
+ 0.5, 0.0, 0.0,
+ 0.0, 1.0, 0.0,
+ 0.5, 1.0, 0.0,
+ 0.0, 0.0, 1.0,
+ 0.5, 0.0, 1.0,
+ 0.0, 1.0, 1.0,
+ 0.5, 1.0, 1.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 1.0, 0.0,
+ 1.0, 0.0, 1.0,
+ 1.0, 1.0, 1.0
+};
+static const size_t nvertices = sizeof(vertices) / sizeof(double[3]);
+
+/* The following array lists the indices toward the 3D vertices of each
+ * triangle.
+ * ,2---,3 ,3---,9 | ,2----3 ,3----9
+ * ,' | ,'/| ,' | ,'/| | ,'/| \ | ,'/| \ | Y
+ * 6----7' / | 7---11' / | | 6' / | \ | 7' / | \ | |
+ * |', | / ,1 |', | / ,8 | | / ,0---,1 | / ,1---,8 o--X
+ * | ',|/,' | ',|/,' | |/,' | ,' |/,' | ,' /
+ * 4----5 5---10 | 4----5' 5---10' Z
+ */
+static const size_t indices[22/*#triangles*/*3/*#indices per triangle*/] = {
+ 0, 4, 2, 2, 4, 6, /* X min */
+ 3, 7, 5, 5, 1, 3, /* X mid */
+ 9,11,10,10, 8, 9, /* X max */
+ 0, 5, 4, 0, 1, 5, 1,10, 5, 1, 8,10, /* Y min */
+ 2, 6, 7, 2, 7, 3, 3, 7,11, 3,11, 9, /* Y max */
+ 0, 2, 1, 1, 2, 3, 1, 3, 8, 8, 3, 9, /* Z min */
+ 4, 5, 6, 6, 5, 7, 5,10, 7, 7,10,11 /* Z max */
+};
+static const size_t ntriangles = sizeof(indices) / sizeof(size_t[3]);
+
+/*******************************************************************************
+ * Box geometry functions
+ ******************************************************************************/
+struct context {
+ const double* vertices;
+ const size_t* indices;
+ struct sdis_interface* interfs[22];
+ const double* scale;
+};
+
+static void
+get_position(const size_t ivert, double pos[3], void* context)
+{
+ struct context* ctx = context;
+ CHK(ctx);
+ pos[0] = ctx->vertices[ivert*3+0] * ctx->scale[0];
+ pos[1] = ctx->vertices[ivert*3+1] * ctx->scale[1];
+ pos[2] = ctx->vertices[ivert*3+2] * ctx->scale[2];
+}
+
+static void
+get_indices(const size_t itri, size_t ids[3], void* context)
+{
+ struct context* ctx = context;
+ CHK(ctx);
+ ids[0] = ctx->indices[itri*3+0];
+ ids[1] = ctx->indices[itri*3+1];
+ ids[2] = ctx->indices[itri*3+2];
+}
+
+static void
+get_interface(const size_t itri, struct sdis_interface** bound, void* context)
+{
+ struct context* ctx = context;
+ CHK(ctx);
+ *bound = ctx->interfs[itri];
+}
+
+/*******************************************************************************
+ * Setup a scene composed of a box that successively contains smaller boxes
+ ******************************************************************************/
+struct matriochka_context {
+ struct sdis_interface* interfs[13];
+ const double* scale;
+ size_t nboxes;
+};
+
+static void
+matriochka_position(const size_t ivert, double pos[3], void* context)
+{
+ struct matriochka_context* ctx = context;
+ const size_t ibox = ctx->nboxes - ivert / box_nvertices - 1;
+ const double* verts = box_vertices;
+ const double box_szmin = 1.0 / (double)ctx->nboxes;
+ const size_t i = ivert % box_nvertices;
+ CHK(ibox <= ctx->nboxes);
+ pos[0] = ((verts[i*3+0]-0.5)*(double)(ibox+1)*box_szmin + 0.5)*ctx->scale[0];
+ pos[1] = ((verts[i*3+1]-0.5)*(double)(ibox+1)*box_szmin + 0.5)*ctx->scale[1];
+ pos[2] = ((verts[i*3+2]-0.5)*(double)(ibox+1)*box_szmin + 0.5)*ctx->scale[2];
+}
+
+static void
+matriochka_indices(const size_t itri, size_t ids[3], void* context)
+{
+ struct matriochka_context* ctx = context;
+ const size_t i = itri % box_ntriangles;
+ const size_t ibox = ctx->nboxes - itri / box_ntriangles - 1;
+ CHK(ibox <= ctx->nboxes);
+ (void)context;
+ ids[0] = box_indices[i*3+0] + ibox*box_nvertices;
+ ids[1] = box_indices[i*3+1] + ibox*box_nvertices;
+ ids[2] = box_indices[i*3+2] + ibox*box_nvertices;
+}
+
+static void
+matriocka_interface
+ (const size_t itri, struct sdis_interface** bound, void* context)
+{
+ struct matriochka_context* ctx = context;
+ const size_t ibox = ctx->nboxes - itri / box_ntriangles - 1;
+ const size_t i = itri % box_ntriangles;
+ CHK(ibox < ctx->nboxes);
+ *bound = ibox != 0 ? ctx->interfs[12] : ctx->interfs[i];
+}
+
+static INLINE void
+dump_matriochkas(FILE* stream, struct matriochka_context* ctx)
+{
+ size_t i;
+ ASSERT(ctx && stream);
+ FOR_EACH(i, 0, ctx->nboxes*box_nvertices) {
+ double pos[3];
+ matriochka_position(i, pos, ctx);
+ fprintf(stream, "v %g %g %g\n", SPLIT3(pos));
+ }
+ FOR_EACH(i, 0, ctx->nboxes*box_ntriangles) {
+ size_t ids[3];
+ matriochka_indices(i, ids, ctx);
+ fprintf(stream, "f %lu %lu %lu\n",
+ (unsigned long)ids[0]+1,
+ (unsigned long)ids[1]+1,
+ (unsigned long)ids[2]+1);
+ }
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double init_temperature;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ if(vtx->time <= 0) {
+ return ((const struct solid*)sdis_data_cget(data))->init_temperature;
+ } else {
+ return UNKNOWN_TEMPERATURE;
+ }
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+struct interf {
+ double temperature;
+};
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->temperature;
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* dev,
+ struct sdis_medium* front,
+ struct sdis_medium* back,
+ const struct sdis_interface_shader* interf_shader,
+ const double temperature)
+{
+ struct sdis_data* data = NULL;
+ struct sdis_interface* interf = NULL;
+ struct interf* interf_props = NULL;
+
+ OK(sdis_data_create
+ (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data));
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = temperature;
+ OK(sdis_interface_create
+ (dev, front, back, interf_shader, data, &interf));
+ OK(sdis_data_ref_put(data));
+ return interf;
+}
+
+/*******************************************************************************
+ * Analytical solution
+ ******************************************************************************/
+static void
+fourier_pq
+ (const size_t nterms_pq,
+ const double pos[3],
+ const double sz[3],
+ const int i0,
+ const int i1,
+ const int i2,
+ double green[2])
+{
+ size_t p, q;
+ CHK(green);
+
+ green[0] = 0;
+ green[1] = 0;
+
+ FOR_EACH(p, 0, nterms_pq+1) {
+ FOR_EACH(q, 0, nterms_pq+1) {
+ const double p2 = (double)(2*p + 1);
+ const double q2 = (double)(2*q + 1);
+ double L_sqr, L, tmp;
+ L_sqr = PI * PI * ((p2*p2)/(sz[i1]*sz[i1]) + (q2*q2)/(sz[i2]*sz[i2]));
+ L = sqrt(L_sqr);
+ tmp = sin(PI*p2*pos[i1]/sz[i1])
+ * sin(PI*q2*pos[i2]/sz[i2])
+ / (sinh(sz[i0]*L)*(p2*q2));
+ if(tmp != 0) {
+ green[0] += sinh(L*(sz[i0]-pos[i0]))*tmp;
+ green[1] += sinh(L*pos[i0])*tmp;
+ }
+ }
+ }
+}
+
+/* This function computes the Green function between a given probe
+ * position/time in a parallelepipedic box and each face of this box (within
+ * the model of a homogeneous boundary condition on each face). */
+static void
+green_analytical
+ (const double box_size[3],
+ const double probe[3],
+ const double time,
+ const double rho,
+ const double cp,
+ const double lambda,
+ double green[7])
+{
+ const size_t nterms_fs = 20; /* #terms in the Fourier expansion series */
+ const size_t nterms_pq = 100; /* #terms in double p/q sums */
+ const size_t nt_pq = (nterms_fs - 1)/2;
+ double Gs[7], Gi[7], Gtmp[7];
+ size_t i, m, n, o, p, q;
+ double a, b, c;
+ double alpha;
+
+ CHK(box_size && probe && time >= 0 && green);
+
+ if(time == 0) {
+ memset(green, 0, sizeof(double[7]));
+ green[6] = 1;
+ return;
+ }
+
+ memset(Gs, 0, sizeof(double[7]));
+ memset(Gi, 0, sizeof(double[7]));
+ memset(Gtmp, 0, sizeof(double[7]));
+
+ /* Steady state solution */
+ fourier_pq(nterms_pq, probe, box_size, 0, 1, 2, Gtmp+0); /* Faces 0 and 1 */
+ fourier_pq(nterms_pq, probe, box_size, 1, 0, 2, Gtmp+2); /* Faces 2 and 3 */
+ fourier_pq(nterms_pq, probe, box_size, 2, 1, 0, Gtmp+4); /* Faces 4 and 5 */
+ FOR_EACH(i, 0, 6) Gs[i] += 16 * Gtmp[i] / (PI * PI);
+
+ alpha=lambda/(rho*cp);
+ a=box_size[0];
+ b=box_size[1];
+ c=box_size[2];
+
+ /* Transient solution */
+ FOR_EACH(m, 0, nterms_fs+1) {
+ const double beta = PI*(double)m/a;
+ const double beta_sqr = beta*beta;
+ const int m_is_even = (m%2 == 0);
+
+ FOR_EACH(n, 0, nterms_fs+1) {
+ const double gamma = PI*(double)n/b;
+ const double gamma_sqr = gamma * gamma;
+ const int n_is_even = (n%2==0);
+
+ FOR_EACH(o, 0, nterms_fs+1) {
+ const double eta = PI*(double)o/c;
+ const double eta_sqr = eta*eta;
+ const double zeta = alpha*(beta_sqr+gamma_sqr+eta_sqr);
+ const int o_is_even = (o%2==0);
+ double Fxyzt;
+
+ memset(Gtmp, 0, sizeof(double[7]));
+ FOR_EACH(p, 0, nt_pq+1) {
+ FOR_EACH(q, 0, nt_pq+1) {
+ const double p2 = (double)(2*p + 1);
+ const double q2 = (double)(2*q + 1);
+ const double Lx_sqr = PI*PI*((p2*p2)/(b*b)+(q2*q2)/(c*c));
+ const double Ly_sqr = PI*PI*((p2*p2)/(a*a)+(q2*q2)/(c*c));
+ const double Lz_sqr = PI*PI*((p2*p2)/(b*b)+(q2*q2)/(a*a));
+ const double pq = p2*q2;
+ double itg[11] = {0};
+
+ itg[1] = 2*q-o+1 == 0 ? -c/2.0 : 0;
+ itg[2] = eta / (eta_sqr+Lz_sqr);
+ itg[4] = 2*p-n+1 == 0 ? -b/2.0 : 0;
+ itg[5] = gamma / (gamma_sqr+Ly_sqr);
+ itg[7] = beta / (beta_sqr+Lx_sqr);
+ itg[9] = 2*p-m+1 == 0 ? -a/2.0 : 0;
+ itg[10] = 2*q-m+1 == 0 ? -a/2.0 : 0;
+
+ if((2*q-o+1==0) && (2*p-n+1==0)) {
+ const double z1 = itg[1]*itg[4]*itg[7];
+ Gtmp[0] += z1/pq;
+ Gtmp[1] += (z1*pow(-1.0,(double)(m+1)))/pq;
+ }
+ if((2*q-o+1==0) && (2*p-m+1==0)) {
+ const double z2 = itg[1]*itg[9]*itg[5];
+ Gtmp[2] += z2/pq;
+ Gtmp[3] += (z2*pow(-1.0,(double)(n+1)))/pq;
+ }
+ if((2*p-n+1==0) && (2*q-m+1==0)) {
+ const double z3 = itg[4]*itg[10]*itg[2];
+ Gtmp[4] += z3/pq;
+ Gtmp[5] += (z3*pow(-1.0,(double)(o+1)))/pq;
+ }
+ }
+ }
+ Fxyzt =
+ sin(probe[0]*beta)
+ * sin(probe[1]*gamma)
+ * sin(probe[2]*eta)
+ * exp(-zeta*time);
+
+ FOR_EACH(i, 0, 6) {
+ Gi[i] += Gtmp[i] * Fxyzt;
+ }
+ if((!m_is_even) && (!n_is_even) && (!o_is_even)) {
+ Gi[6] += 8.0 * Fxyzt/(beta*gamma*eta);
+ }
+ }
+ }
+ }
+
+ /* Gi[i], i=0,5: Green of boundary index i */
+ FOR_EACH(i, 0, 6) {
+ Gi[i] = -128.0 * Gi[i]/(a*b*c*PI*PI);
+ }
+
+ /* Gi[6]: Green of the initial condition */
+ Gi[6] = 8.0*Gi[6]/(a*b*c);
+
+ /* Computing total Green function */
+ FOR_EACH(i, 0, 7) {
+ green[i] = Gs[i] + Gi[i];
+ }
+}
+
+static double
+temperature_analytical
+ (const double temperature_bounds[6],
+ const double temperature_init,
+ const double box_size[3],
+ const double probe[3],
+ const double time,
+ const double rho,
+ const double cp,
+ const double lambda)
+{
+ double green[7];
+ double temperature = 0;
+ size_t i;
+ CHK(temperature_bounds && temperature_init && box_size && probe);
+ green_analytical(box_size, probe, time, rho, cp, lambda, green);
+
+ FOR_EACH(i, 0, 6) {
+ printf("Green for face %lu: %g\n", (unsigned long)i, green[i]);
+ temperature += green[i] * temperature_bounds[i];
+ }
+ temperature += green[6] * temperature_init;
+ return temperature;
+}
+
+/*******************************************************************************
+ * Main function
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_device* dev = NULL;
+ struct sdis_scene* box_scn = NULL;
+ struct sdis_scene* box2_scn = NULL;
+ struct sdis_scene* box_matriochka_scn = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* interfs[7] = {NULL};
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_mc temperature = SDIS_MC_NULL;
+ struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct solid* solid_param = NULL;
+ struct context ctx;
+ struct matriochka_context matriochka_ctx;
+ const size_t nrealisations = 10000;
+ const size_t nmatriochkas = 5;
+ size_t nfails = 0;
+ double probe[3];
+ double time[2];
+ double Tbounds[6];
+ double Tinit;
+ double Tref;
+ double boxsz[3];
+ double rho, cp, lambda;
+ (void)argc, (void)argv;
+
+ /* System description */
+ rho = 1700;
+ cp = 800;
+ lambda = 1.15;
+ Tbounds[0] = 280; /* Xmin */
+ Tbounds[1] = 290; /* Xmax */
+ Tbounds[2] = 310; /* Ymin */
+ Tbounds[3] = 270; /* Ymax */
+ Tbounds[4] = 300; /* Zmin */
+ Tbounds[5] = 320; /* Zmax */
+ Tinit = 100;
+ boxsz[0] = 0.3;
+ boxsz[1] = 0.1;
+ boxsz[2] = 0.2;
+
+ OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+
+ /* Create the fluid medium */
+ fluid_shader = DUMMY_FLUID_SHADER;
+ OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
+
+ /* 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 = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+
+ /* 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->rho = rho;
+ solid_param->cp = cp;
+ solid_param->lambda = lambda;
+ solid_param->delta = 1.0/20.0 * MMIN(MMIN(boxsz[0], boxsz[1]), boxsz[2]);
+ solid_param->init_temperature = Tinit;
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid));
+
+ /* Setup the interface shader */
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the interfaces */
+ interfs[0] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[0]);
+ interfs[1] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[1]);
+ interfs[2] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[2]);
+ interfs[3] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[3]);
+ interfs[4] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[4]);
+ interfs[5] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[5]);
+ interfs[6] = create_interface(dev, solid, solid, &interf_shader, UNKNOWN_TEMPERATURE);
+
+ /* Setup the box scene context */
+ ctx.indices = box_indices;
+ ctx.vertices = box_vertices;
+ ctx.interfs[0] = ctx.interfs[1] = interfs[4]; /* Zmin */
+ ctx.interfs[2] = ctx.interfs[3] = interfs[0]; /* Xmin */
+ ctx.interfs[4] = ctx.interfs[5] = interfs[5]; /* Zmax */
+ ctx.interfs[6] = ctx.interfs[7] = interfs[1]; /* Xmax */
+ ctx.interfs[8] = ctx.interfs[9] = interfs[3]; /* Ymax */
+ ctx.interfs[10] = ctx.interfs[11] = interfs[2]; /* Ymin */
+ ctx.scale = boxsz;
+
+ /* Create the box scene */
+ OK(sdis_scene_create(dev, box_ntriangles, get_indices, get_interface,
+ box_nvertices, get_position, &ctx, &box_scn));
+
+ /* Setup the box2 scene context */
+ ctx.indices = indices;
+ ctx.vertices = vertices;
+ ctx.interfs[0] = ctx.interfs[1] = interfs[0]; /* Xmin */
+ ctx.interfs[2] = ctx.interfs[3] = interfs[6]; /* Xmid */
+ ctx.interfs[4] = ctx.interfs[5] = interfs[1]; /* Xmax */
+ ctx.interfs[6] = ctx.interfs[7] = interfs[2]; /* Ymin */
+ ctx.interfs[8] = ctx.interfs[9] = interfs[2]; /* Ymin */
+ ctx.interfs[10] = ctx.interfs[11] = interfs[3]; /* Ymax */
+ ctx.interfs[12] = ctx.interfs[13] = interfs[3]; /* Ymax */
+ ctx.interfs[14] = ctx.interfs[15] = interfs[4]; /* Zmin */
+ ctx.interfs[16] = ctx.interfs[17] = interfs[4]; /* Zmin */
+ ctx.interfs[18] = ctx.interfs[19] = interfs[5]; /* Zmax */
+ ctx.interfs[20] = ctx.interfs[21] = interfs[5]; /* Zmax */
+ ctx.scale = boxsz;
+
+ /* Create the box scene */
+ OK(sdis_scene_create(dev, ntriangles, get_indices, get_interface,
+ nvertices, get_position, &ctx, &box2_scn));
+
+ /* Setup the matriochka context */
+ matriochka_ctx.interfs[0] = matriochka_ctx.interfs[1] = interfs[4]; /* Zmin */
+ matriochka_ctx.interfs[2] = matriochka_ctx.interfs[3] = interfs[0]; /* Xmin */
+ matriochka_ctx.interfs[4] = matriochka_ctx.interfs[5] = interfs[5]; /* Zmax */
+ matriochka_ctx.interfs[6] = matriochka_ctx.interfs[7] = interfs[1]; /* Xmax */
+ matriochka_ctx.interfs[8] = matriochka_ctx.interfs[9] = interfs[3]; /* Ymax */
+ matriochka_ctx.interfs[10] = matriochka_ctx.interfs[11] = interfs[2]; /* Ymin */
+ matriochka_ctx.interfs[12] = interfs[6]; /* The remaining internal triangles */
+ matriochka_ctx.scale = boxsz;
+ matriochka_ctx.nboxes = nmatriochkas;
+
+ /* Create the matriochka scene */
+ OK(sdis_scene_create(dev, box_ntriangles*nmatriochkas, matriochka_indices,
+ matriocka_interface, box_nvertices*nmatriochkas, matriochka_position,
+ &matriochka_ctx, &box_matriochka_scn));
+
+ /* Setup and run the simulation */
+ probe[0] = 0.1;
+ probe[1] = 0.06;
+ probe[2] = 0.130;
+ time[0] = time[1] = 500; /* Observation time range */
+
+ /* Compute the solution */
+ Tref = temperature_analytical
+ (Tbounds, Tinit, boxsz, probe, time[0], rho, cp, lambda);
+
+ /* Run simulation on regular scene */
+ solid_param->delta = 1.0/20.0 * MMIN(MMIN(boxsz[0], boxsz[1]), boxsz[2]);
+ solve_args.nrealisations = nrealisations;
+ solve_args.position[0] = probe[0];
+ solve_args.position[1] = probe[1];
+ solve_args.position[2] = probe[2];
+ solve_args.time_range[0] = time[0];
+ solve_args.time_range[1] = time[1];
+ OK(sdis_solve_probe(box_scn, &solve_args, &estimator));
+
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &temperature));
+ printf("Temperature at (%g, %g, %g) m at %g s (delta = %g) = %g ~ %g +/- %g\n",
+ SPLIT3(probe), time[0], solid_param->delta, Tref, temperature.E,
+ temperature.SE);
+ printf("#failures = %lu/%lu\n",
+ (unsigned long)nfails, (unsigned long)nrealisations);
+ CHK(eq_eps(Tref, temperature.E, temperature.SE*3));
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Run simulation on split scene */
+ solid_param->delta = 1.0/20.0 * MMIN(MMIN(boxsz[0]/2.0, boxsz[1]), boxsz[2]);
+ OK(sdis_solve_probe(box2_scn, &solve_args, &estimator));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &temperature));
+ printf("Temperature at (%g, %g, %g) m at %g s (delta = %g) = %g ~ %g +/- %g\n",
+ SPLIT3(probe), time[0], solid_param->delta, Tref, temperature.E,
+ temperature.SE);
+ printf("#failures = %lu/%lu\n",
+ (unsigned long)nfails, (unsigned long)nrealisations);
+ CHK(eq_eps(Tref, temperature.E, temperature.SE*3));
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Run simulation on matriochkas */
+ solid_param->delta = MMIN(MMIN(boxsz[0], boxsz[1]), boxsz[2]);
+ solid_param->delta /= (double)nmatriochkas;
+ solid_param->delta *= 1.0/20.0;
+ OK(sdis_solve_probe(box_matriochka_scn, &solve_args, &estimator));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &temperature));
+ printf("Temperature at (%g, %g, %g) m at %g s (delta = %g) = %g ~ %g +/- %g\n",
+ SPLIT3(probe), time[0], solid_param->delta, Tref, temperature.E,
+ temperature.SE);
+ printf("#failures = %lu/%lu\n",
+ (unsigned long)nfails, (unsigned long)nrealisations);
+ CHK(eq_eps(Tref, temperature.E, temperature.SE*3));
+ OK(sdis_estimator_ref_put(estimator));
+
+ OK(sdis_interface_ref_put(interfs[0]));
+ OK(sdis_interface_ref_put(interfs[1]));
+ OK(sdis_interface_ref_put(interfs[2]));
+ OK(sdis_interface_ref_put(interfs[3]));
+ OK(sdis_interface_ref_put(interfs[4]));
+ OK(sdis_interface_ref_put(interfs[5]));
+ OK(sdis_interface_ref_put(interfs[6]));
+ OK(sdis_data_ref_put(data));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_medium_ref_put(fluid));
+ OK(sdis_device_ref_put(dev));
+ OK(sdis_scene_ref_put(box_scn));
+ OK(sdis_scene_ref_put(box2_scn));
+ OK(sdis_scene_ref_put(box_matriochka_scn));
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -363,3 +363,35 @@ dump_heat_paths(FILE* stream, const struct sdis_estimator* estimator)
fprintf(stream, "0.0 0.0 1.0 1.0\n"); /* 0.0 = Blue: success */
fprintf(stream, "1.0 0.0 0.0 1.0\n"); /* 1.0 = Red: failure */
}
+
+void
+check_green_serialization
+ (struct sdis_green_function* green,
+ struct sdis_scene* scn,
+ const double time_range[2])
+{
+ FILE* stream = NULL;
+ struct sdis_estimator *e1 = NULL;
+ struct sdis_estimator *e2 = NULL;
+ struct sdis_green_function* green2 = NULL;
+
+ CHK(green && time_range);
+ CHK(stream = tmpfile());
+
+ OK(sdis_green_function_write(green, stream));
+
+ rewind(stream);
+ OK(sdis_green_function_create_from_stream(scn, stream, &green2));
+ CHK(!fclose(stream));
+ check_green_function(green2);
+
+ OK(sdis_green_function_solve(green, time_range, &e1));
+ OK(sdis_green_function_solve(green2, time_range, &e2));
+ check_estimator_eq_strict(e1, e2);
+
+ OK(sdis_estimator_ref_put(e1));
+ OK(sdis_estimator_ref_put(e2));
+ OK(sdis_green_function_ref_put(green2));
+}
+
+
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -18,6 +18,7 @@
#include "sdis.h"
+#include <rsys/double33.h>
#include <rsys/mem_allocator.h>
#include <stdio.h>
@@ -52,12 +53,12 @@ static const size_t box_nvertices = sizeof(box_vertices) / sizeof(double[3]);
* Front, right Back, left and Z
* and Top faces bottom faces */
static const size_t box_indices[12/*#triangles*/*3/*#indices per triangle*/] = {
- 0, 2, 1, 1, 2, 3, /* Front face */
- 0, 4, 2, 2, 4, 6, /* Left face*/
- 4, 5, 6, 6, 5, 7, /* Back face */
- 3, 7, 1, 1, 7, 5, /* Right face */
- 2, 6, 3, 3, 6, 7, /* Top face */
- 0, 1, 4, 4, 1, 5 /* Bottom face */
+ 0, 2, 1, 1, 2, 3, /* -Z */
+ 0, 4, 2, 2, 4, 6, /* -X */
+ 4, 5, 6, 6, 5, 7, /* +Z */
+ 3, 7, 5, 5, 1, 3, /* +X */
+ 2, 6, 7, 7, 3, 2, /* +Y */
+ 0, 1, 5, 5, 4, 0 /* -Y */
};
static const size_t box_ntriangles = sizeof(box_indices) / sizeof(size_t[3]);
@@ -251,26 +252,82 @@ check_estimator_eq
OK(sdis_estimator_get_type(e2, &type2));
CHK(type1 == type2);
- OK(sdis_estimator_get_temperature(e1, &mc1));
- OK(sdis_estimator_get_temperature(e2, &mc2));
- CHK(mc1.E + 3*mc1.SE >= mc2.E - 3*mc2.SE);
- CHK(mc1.E - 3*mc1.SE <= mc2.E + 3*mc2.SE);
-
- if(type1 == SDIS_ESTIMATOR_FLUX) {
- OK(sdis_estimator_get_convective_flux(e1, &mc1));
- OK(sdis_estimator_get_convective_flux(e2, &mc2));
- CHK(mc1.E + 3*mc1.SE >= mc2.E - 3*mc2.SE);
- CHK(mc1.E - 3*mc1.SE <= mc2.E + 3*mc2.SE);
-
- OK(sdis_estimator_get_radiative_flux(e1, &mc1));
- OK(sdis_estimator_get_radiative_flux(e2, &mc2));
- CHK(mc1.E + 3*mc1.SE >= mc2.E - 3*mc2.SE);
- CHK(mc1.E - 3*mc1.SE <= mc2.E + 3*mc2.SE);
-
- OK(sdis_estimator_get_total_flux(e1, &mc1));
- OK(sdis_estimator_get_total_flux(e2, &mc2));
- CHK(mc1.E + 3*mc1.SE >= mc2.E - 3*mc2.SE);
- CHK(mc1.E - 3*mc1.SE <= mc2.E + 3*mc2.SE);
+ switch(type1) {
+ case SDIS_ESTIMATOR_TEMPERATURE:
+ OK(sdis_estimator_get_temperature(e1, &mc1));
+ OK(sdis_estimator_get_temperature(e2, &mc2));
+ CHK(mc1.E + 3*mc1.SE >= mc2.E - 3*mc2.SE);
+ CHK(mc1.E - 3*mc1.SE <= mc2.E + 3*mc2.SE);
+ break;
+
+ case SDIS_ESTIMATOR_FLUX:
+ OK(sdis_estimator_get_convective_flux(e1, &mc1));
+ OK(sdis_estimator_get_convective_flux(e2, &mc2));
+ CHK(mc1.E + 3*mc1.SE >= mc2.E - 3*mc2.SE);
+ CHK(mc1.E - 3*mc1.SE <= mc2.E + 3*mc2.SE);
+
+ OK(sdis_estimator_get_radiative_flux(e1, &mc1));
+ OK(sdis_estimator_get_radiative_flux(e2, &mc2));
+ CHK(mc1.E + 3*mc1.SE >= mc2.E - 3*mc2.SE);
+ CHK(mc1.E - 3*mc1.SE <= mc2.E + 3*mc2.SE);
+
+ OK(sdis_estimator_get_total_flux(e1, &mc1));
+ OK(sdis_estimator_get_total_flux(e2, &mc2));
+ CHK(mc1.E + 3*mc1.SE >= mc2.E - 3*mc2.SE);
+ CHK(mc1.E - 3*mc1.SE <= mc2.E + 3*mc2.SE);
+ break;
+
+ case SDIS_ESTIMATOR_POWER:
+ OK(sdis_estimator_get_power(e1, &mc1));
+ OK(sdis_estimator_get_power(e2, &mc2));
+ CHK(mc1.E + 3*mc1.SE >= mc2.E - 3*mc2.SE);
+ CHK(mc1.E - 3*mc1.SE <= mc2.E + 3*mc2.SE);
+ break;
+
+ default: FATAL("Unreachable code.\n"); break;
+ }
+}
+
+static INLINE void
+check_estimator_eq_strict
+ (const struct sdis_estimator* e1, const struct sdis_estimator* e2)
+{
+ struct sdis_mc mc1, mc2;
+ enum sdis_estimator_type type1, type2;
+ ASSERT(e1 && e2);
+
+ OK(sdis_estimator_get_type(e1, &type1));
+ OK(sdis_estimator_get_type(e2, &type2));
+ CHK(type1 == type2);
+
+ switch(type1) {
+ case SDIS_ESTIMATOR_TEMPERATURE:
+ OK(sdis_estimator_get_temperature(e1, &mc1));
+ OK(sdis_estimator_get_temperature(e2, &mc2));
+ CHK(mc1.E == mc2.E && mc1.V == mc2.V && mc1.SE == mc2.SE);
+ break;
+
+ case SDIS_ESTIMATOR_FLUX:
+ OK(sdis_estimator_get_convective_flux(e1, &mc1));
+ OK(sdis_estimator_get_convective_flux(e2, &mc2));
+ CHK(mc1.E == mc2.E && mc1.V == mc2.V && mc1.SE == mc2.SE);
+
+ OK(sdis_estimator_get_radiative_flux(e1, &mc1));
+ OK(sdis_estimator_get_radiative_flux(e2, &mc2));
+ CHK(mc1.E == mc2.E && mc1.V == mc2.V && mc1.SE == mc2.SE);
+
+ OK(sdis_estimator_get_total_flux(e1, &mc1));
+ OK(sdis_estimator_get_total_flux(e2, &mc2));
+ CHK(mc1.E == mc2.E && mc1.V == mc2.V && mc1.SE == mc2.SE);
+ break;
+
+ case SDIS_ESTIMATOR_POWER:
+ OK(sdis_estimator_get_power(e1, &mc1));
+ OK(sdis_estimator_get_power(e2, &mc2));
+ CHK(mc1.E == mc2.E && mc1.V == mc2.V && mc1.SE == mc2.SE);
+ break;
+
+ default: FATAL("Unreachable code.\n"); break;
}
}
@@ -294,5 +351,11 @@ dump_heat_paths
(FILE* stream,
const struct sdis_estimator* estimator);
+extern LOCAL_SYM void
+check_green_serialization
+ (struct sdis_green_function* green,
+ struct sdis_scene* scn,
+ const double time_range[2]);
+
#endif /* TEST_SDIS_UTILS_H */
diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c
@@ -60,6 +60,7 @@ struct solid {
double rho;
double cp;
double delta;
+ double vpower;
};
static double
@@ -118,7 +119,7 @@ solid_get_volumic_power
{
(void)data;
CHK(vtx != NULL);
- return P0;
+ return ((struct solid*)sdis_data_cget(data))->vpower;
}
/*******************************************************************************
@@ -149,7 +150,10 @@ interface_get_convection_coef
* Helper functions
******************************************************************************/
static void
-solve(struct sdis_scene* scn, const double pos[])
+solve
+ (struct sdis_scene* scn,
+ const double pos[],
+ struct solid* solid)
{
char dump[128];
struct time t0, t1, t2;
@@ -157,17 +161,21 @@ solve(struct sdis_scene* scn, const double pos[])
struct sdis_estimator* estimator2;
struct sdis_green_function* green;
struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
- struct sdis_mc T;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
size_t nreals;
size_t nfails;
double x;
double ref;
const double time_range[2] = {INF, INF};
enum sdis_scene_dimension dim;
- ASSERT(scn && pos);
+ ASSERT(scn && pos && solid);
+
+ /* Restore power value */
+ solid->vpower = P0;
x = pos[0] - 0.5;
- ref = P0 / (2*LAMBDA) * (1.0/4.0 - x*x) + T0;
+ ref = solid->vpower / (2*LAMBDA) * (1.0/4.0 - x*x) + T0;
OK(sdis_scene_get_dimension(scn, &dim));
@@ -186,25 +194,28 @@ solve(struct sdis_scene* scn, const double pos[])
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("Temperature at (%g %g) = %g ~ %g +/- %g [%g, %g]\n",
- SPLIT2(pos), ref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE);
+ printf("Temperature at (%g %g) with Power=%g = %g ~ %g +/- %g [%g, %g]\n",
+ SPLIT2(pos), solid->vpower, ref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE);
break;
case SDIS_SCENE_3D:
- printf("Temperature at (%g %g %g) = %g ~ %g +/- %g [%g, %g]\n",
- SPLIT3(pos), ref, T.E, T.SE, T.E -3*T.SE, T.E + 3*T.SE);
+ printf("Temperature at (%g %g %g)th Power=%g = %g ~ %g +/- %g [%g, %g]\n",
+ SPLIT3(pos), solid->vpower, ref, T.E, T.SE, T.E -3*T.SE, T.E + 3*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\n", dump);
+ 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, T.SE*3));
+ /* Check green function */
time_current(&t0);
OK(sdis_solve_probe_green_function(scn, &solve_args, &green));
time_current(&t1);
@@ -217,12 +228,12 @@ solve(struct sdis_scene* scn, const double pos[])
switch(dim) {
case SDIS_SCENE_2D:
- printf("Green temperature at (%g %g) = %g ~ %g +/- %g\n",
- SPLIT2(pos), ref, T.E, T.SE);
+ printf("Green temperature at (%g %g) with Power=%g = %g ~ %g +/- %g\n",
+ SPLIT2(pos), solid->vpower, ref, T.E, T.SE);
break;
case SDIS_SCENE_3D:
- printf("Green temperature at (%g %g %g) = %g ~ %g +/- %g\n",
- SPLIT3(pos), ref, T.E, T.SE);
+ printf("Green temperature at (%g %g %g) with Power=%g = %g ~ %g +/- %g\n",
+ SPLIT3(pos), solid->vpower, ref, T.E, T.SE);
break;
default: FATAL("Unreachable code.\n"); break;
}
@@ -236,6 +247,54 @@ solve(struct sdis_scene* scn, const double pos[])
check_green_function(green);
check_estimator_eq(estimator, estimator2);
+ check_green_serialization(green, scn, time_range);
+
+ OK(sdis_estimator_ref_put(estimator));
+ OK(sdis_estimator_ref_put(estimator2));
+ printf("\n");
+
+ /* Check green used at a different power */
+ solid->vpower = 3 * P0;
+
+ 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));
+
+ ref = solid->vpower / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0;
+
+ switch (dim) {
+ case SDIS_SCENE_2D:
+ printf("Temperature at (%g %g) with Power=%g = %g ~ %g +/- %g [%g, %g]\n",
+ SPLIT2(pos), solid->vpower, ref, T.E, T.SE, T.E - 3 * T.SE, T.E + 3 * T.SE);
+ break;
+ case SDIS_SCENE_3D:
+ printf("Temperature at (%g %g %g) with Power=%g = %g ~ %g +/- %g [%g, %g]\n",
+ SPLIT3(pos), solid->vpower, ref, T.E, T.SE, T.E - 3 * T.SE, T.E + 3 * 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", time.E, time.SE);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails < N / 1000);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+
+ time_current(&t0);
+ OK(sdis_green_function_solve(green, time_range, &estimator2));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+ printf("Green solve time = %s\n", dump);
+
+ check_green_function(green);
+ check_estimator_eq(estimator, estimator2);
OK(sdis_estimator_ref_put(estimator));
OK(sdis_estimator_ref_put(estimator2));
@@ -255,7 +314,6 @@ main(int argc, char** argv)
struct sdis_device* dev = NULL;
struct sdis_medium* fluid = NULL;
struct sdis_medium* solid = NULL;
- struct sdis_medium* solid2 = NULL; /* For debug */
struct sdis_interface* interf_adiabatic = NULL;
struct sdis_interface* interf_T0 = NULL;
struct sdis_scene* box_scn = NULL;
@@ -291,18 +349,10 @@ main(int argc, char** argv)
solid_props->cp = 2;
solid_props->rho = 25;
solid_props->delta = DELTA;
+ solid_props->vpower = P0;
OK(sdis_solid_create(dev, &solid_shader, data, &solid));
OK(sdis_data_ref_put(data));
- OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data));
- solid_props = sdis_data_get(data);
- solid_props->lambda = 0;
- solid_props->cp = 0;
- solid_props->rho = 0;
- solid_props->delta = DELTA/4;
- OK(sdis_solid_create(dev, &solid_shader, data, &solid2));
- OK(sdis_data_ref_put(data));
-
/* Setup the interface shader */
interf_shader.convection_coef = interface_get_convection_coef;
interf_shader.front.temperature = interface_get_temperature;
@@ -325,7 +375,6 @@ main(int argc, char** argv)
/* Release the media */
OK(sdis_medium_ref_put(solid));
- OK(sdis_medium_ref_put(solid2));
OK(sdis_medium_ref_put(fluid));
/* Map the interfaces to their box triangles */
@@ -359,9 +408,9 @@ main(int argc, char** argv)
/* Solve */
d3_splat(pos, 0.25);
printf(">> Box scene\n");
- solve(box_scn, pos);
+ solve(box_scn, pos, solid_props);
printf(">> Square scene\n");
- solve(square_scn, pos);
+ solve(square_scn, pos, solid_props);
OK(sdis_scene_ref_put(box_scn));
OK(sdis_scene_ref_put(square_scn));
