commit 1a2b6ef0f7459e2541e56be74ef2a7749fd1d52d
parent 3fc7ccb62e9be7ae9821810680fa6fbdcdbba9db
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 17 Jul 2020 15:57:55 +0200
Merge branch 'feature_mean_power' into develop
Diffstat:
7 files changed, 621 insertions(+), 5 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -172,6 +172,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)
@@ -208,6 +209,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
@@ -113,8 +113,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_MEAN_POWER, /* In Watt */
SDIS_ESTIMATOR_TYPES_COUNT__
};
@@ -501,6 +502,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_mean_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_MEAN_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_mean_power_args
+SDIS_COMPUTE_MEAN_POWER_ARGS_DEFAULT = SDIS_COMPUTE_MEAN_POWER_ARGS_DEFAULT__;
+
BEGIN_DECLS
/*******************************************************************************
@@ -912,6 +930,11 @@ sdis_estimator_get_total_flux
struct sdis_mc* flux);
SDIS_API res_T
+sdis_estimator_get_mean_power
+ (const struct sdis_estimator* estimator,
+ struct sdis_mc* mean_power);
+
+SDIS_API res_T
sdis_estimator_get_paths_count
(const struct sdis_estimator* estimator,
size_t* npaths);
@@ -1083,6 +1106,15 @@ sdis_solve_medium
const struct sdis_solve_medium_args* args,
struct sdis_estimator** estimator);
+/* P = SUM(volumic_power(x)) / Nrealisations * Volume
+ * mean power (in Watt) = time_range[0] == time_range[1]
+ * ? P : P / (time_range[1] - time_range[0]) */
+SDIS_API res_T
+sdis_compute_mean_power
+ (struct sdis_scene* scn,
+ const struct sdis_compute_mean_power_args* args,
+ struct sdis_estimator** estimator);
+
/*******************************************************************************
* Green solvers.
*
diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c
@@ -118,7 +118,6 @@ sdis_estimator_get_convective_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 +128,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 +138,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_mean_power
+ (const struct sdis_estimator* estimator, struct sdis_mc* mean_power)
+{
+ if(!estimator || !mean_power || estimator->type != SDIS_ESTIMATOR_MEAN_POWER)
+ return RES_BAD_ARG;
+ SETUP_MC(mean_power, &estimator->mean_power.power);
+ mean_power->E *= estimator->mean_power.spread;
+
+ if(estimator->mean_power.time_range[0]
+ != estimator->mean_power.time_range[1]) {
+ mean_power->E /= /* From Joule to Watt */
+ (estimator->mean_power.time_range[1]-estimator->mean_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];
+ } mean_power;
+
size_t nrealisations; /* #successes */
size_t nfailures;
@@ -98,6 +105,23 @@ estimator_setup_temperature
}
static INLINE void
+estimator_setup_mean_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->mean_power.power.sum = sum;
+ estim->mean_power.power.sum2 = sum2;
+ estim->mean_power.power.count = estim->nrealisations;
+ estim->mean_power.spread = spread;
+ estim->mean_power.time_range[0] = time_range[0];
+ estim->mean_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_solve.c b/src/sdis_solve.c
@@ -540,3 +540,16 @@ sdis_solve_medium_green_function
}
}
+res_T
+sdis_compute_mean_power
+ (struct sdis_scene* scn,
+ const struct sdis_compute_mean_power_args* args,
+ struct sdis_estimator** estimator)
+{
+ if(!scn) return RES_BAD_ARG;
+ if(scene_is_2d(scn)) {
+ return compute_mean_power_2d(scn, args, estimator);
+ } else {
+ return compute_mean_power_3d(scn, args, estimator);
+ }
+}
diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h
@@ -477,4 +477,192 @@ error:
goto exit;
}
+static res_T
+XD(compute_mean_power)
+ (struct sdis_scene* scn,
+ const struct sdis_compute_mean_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_MEAN_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_mean_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/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_mean_power_args args = SDIS_COMPUTE_MEAN_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_mean_power function */
+ args.nrealisations = N;
+ args.medium = solid0;
+ args.time_range[0] = INF;
+ args.time_range[1] = INF;
+ BA(sdis_compute_mean_power(NULL, &args, &estimator));
+ BA(sdis_compute_mean_power(scn, NULL, &estimator));
+ BA(sdis_compute_mean_power(scn, &args, NULL));
+ args.nrealisations = 0;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+ args.nrealisations = N;
+ args.medium = NULL;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+ args.medium = solid0;
+ args.fp_to_meter = 0;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+ args.fp_to_meter = 1;
+ args.time_range[0] = args.time_range[1] = -1;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+ args.time_range[0] = 1;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+ args.time_range[1] = 0;
+ BA(sdis_compute_mean_power(scn, &args, &estimator));
+ args.time_range[0] = args.time_range[1] = INF;
+ OK(sdis_compute_mean_power(scn, &args, &estimator));
+
+ BA(sdis_estimator_get_mean_power(NULL, &mpow));
+ BA(sdis_estimator_get_mean_power(estimator, NULL));
+ OK(sdis_estimator_get_mean_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_mean_power(scn, &args, &estimator));
+ OK(sdis_estimator_get_mean_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_mean_power(scn, &args, &estimator));
+ OK(sdis_estimator_get_mean_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_mean_power(scn, &args, &estimator));
+
+ /* Check non constant volumic power */
+ args.medium = solid0;
+ OK(sdis_compute_mean_power(scn, &args, &estimator));
+ OK(sdis_estimator_get_mean_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.e-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;
+}
