commit 37b228efacad898be2bcfc4f3a9c63910516208d
parent 03c8f1052ed31a29fbd48725f60019db488773ef
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Tue, 17 Jan 2017 18:30:01 +0100
First version of global MC result by receiver.
Diffstat:
10 files changed, 184 insertions(+), 15 deletions(-)
diff --git a/src/ssol.h b/src/ssol.h
@@ -754,6 +754,13 @@ ssol_estimator_get_status
struct ssol_estimator_status* status);
SSOL_API res_T
+ssol_estimator_get_receiver_status
+ (const struct ssol_estimator* estimator,
+ const struct ssol_instance* instance,
+ const enum ssol_side_flag side,
+ struct ssol_estimator_status* status);
+
+SSOL_API res_T
ssol_estimator_get_count
(const struct ssol_estimator* estimator,
size_t* count);
diff --git a/src/ssol_estimator.c b/src/ssol_estimator.c
@@ -16,6 +16,7 @@
#include "ssol.h"
#include "ssol_estimator_c.h"
#include "ssol_device_c.h"
+#include "ssol_instance_c.h"
#include <rsys/rsys.h>
#include <rsys/mem_allocator.h>
@@ -34,6 +35,7 @@ estimator_release(ref_T* ref)
CONTAINER_OF(ref, struct ssol_estimator, ref);
ASSERT(ref);
dev = estimator->dev;
+ darray_receiver_release(&estimator->global_receiver);
ASSERT(dev && dev->allocator);
MEM_RM(dev->allocator, estimator);
SSOL(device_ref_put(dev));
@@ -61,6 +63,8 @@ ssol_estimator_create
goto error;
}
+ darray_receiver_init(dev->allocator, &estimator->global_receiver);
+
SSOL(device_ref_get(dev));
estimator->dev = dev;
ref_init(&estimator->ref);
@@ -94,6 +98,7 @@ ssol_estimator_ref_put
ref_put(&estimator->ref, estimator_release);
return RES_OK;
}
+
res_T
ssol_estimator_get_status
(const struct ssol_estimator* estimator,
@@ -118,6 +123,34 @@ ssol_estimator_get_status
}
res_T
+ssol_estimator_get_receiver_status
+ (const struct ssol_estimator* estimator,
+ const struct ssol_instance* instance,
+ const enum ssol_side_flag side,
+ struct ssol_estimator_status* status)
+{
+ const struct mc_data* data;
+ uint32_t idx = 0;
+ if (!estimator || !instance || !status
+ || (side != SSOL_BACK && side != SSOL_FRONT))
+ return RES_BAD_ARG;
+
+ /* check if a receiver is defined for this instance/side */
+ idx = instance->mc_result_idx[side_idx(side)];
+ if (idx == MC_RESULT_IDX_NONE)
+ return RES_BAD_ARG;
+
+ ASSERT(idx < estimator->global_receiver.size);
+ data = estimator->global_receiver.data + idx;
+ status->N = estimator->realisation_count;
+ status->Nf = estimator->failed_count;
+ status->E = data->weight / (double) status->N;
+ status->V = data->sqr_weight / (double) status->N - status->E * status->E;
+ status->SE = (status->V > 0) ? sqrt(status->V / (double) status->N) : 0;
+ return RES_OK;
+}
+
+res_T
ssol_estimator_get_count
(const struct ssol_estimator* estimator,
size_t* count)
diff --git a/src/ssol_estimator_c.h b/src/ssol_estimator_c.h
@@ -17,6 +17,9 @@
#define SSOL_ESTIMATOR_C_H
#include <rsys/ref_count.h>
+#include <rsys/dynamic_array.h>
+
+struct mem_allocator;
/* Monte carlo data */
struct mc_data {
@@ -29,13 +32,30 @@ struct mc_data {
#define MC_DATA_NULL__ { 0, 0 }
static const struct mc_data MC_DATA_NULL = MC_DATA_NULL__;
+static INLINE void
+init_mc_data
+ (struct mem_allocator* alloc,
+ struct mc_data* data)
+{
+ (void)alloc;
+ ASSERT(data);
+ CLEAR_MC_DATA(*data);
+}
+
+/* Define the darray_receiver data structure */
+#define DARRAY_NAME receiver
+#define DARRAY_DATA struct mc_data
+#define DARRAY_FUNCTOR_INIT init_mc_data
+#include <rsys/dynamic_array.h>
+
struct ssol_estimator {
size_t realisation_count;
size_t failed_count;
/* the implicit MC computations */
struct mc_data shadow;
struct mc_data missing;
- /* 2 global MC per receiver: one for P, one for cos effect losses */
+ /* 1 global MC per receiver */
+ struct darray_receiver global_receiver;
struct ssol_device* dev;
ref_T ref;
diff --git a/src/ssol_instance.c b/src/ssol_instance.c
@@ -78,6 +78,8 @@ ssol_object_instantiate
instance->dev = dev;
instance->object = object;
instance->sample = 1;
+ instance->mc_result_idx[side_idx(SSOL_FRONT)] = MC_RESULT_IDX_NONE;
+ instance->mc_result_idx[side_idx(SSOL_BACK)] = MC_RESULT_IDX_NONE;
d33_set_identity(instance->transform);
d3_splat(instance->transform + 9, 0);
diff --git a/src/ssol_instance_c.h b/src/ssol_instance_c.h
@@ -20,13 +20,23 @@
#include <rsys/list.h>
#include <rsys/ref_count.h>
+#define MC_RESULT_IDX_NONE -1
+
+static FINLINE int
+side_idx(const enum ssol_side_flag side)
+{
+ ASSERT(side == SSOL_FRONT || side == SSOL_BACK);
+ return side == SSOL_FRONT ? 0 : 1;
+}
+
struct ssol_instance {
struct ssol_object* object; /* Instantiated object */
struct s3d_shape* shape_rt; /* Instantiated Star-3D shape to ray-trace */
struct s3d_shape* shape_samp; /* Instantiated Star-3D shape to sample */
double transform[12]; /* Column major 4x3 affine transformation */
+ int mc_result_idx[2]; /* index of the global MC results for possible receptors */
int receiver_mask; /* Combination of ssol_face_flag */
- int sample; /* Define whether or not the instance is sampled or not */
+ int sample; /* Define whether or not the instance should be sampled */
struct fid id; /* Unique identifier */
diff --git a/src/ssol_scene.c b/src/ssol_scene.c
@@ -282,7 +282,7 @@ scene_create_s3d_views
struct s3d_scene_view* view_rt = NULL;
struct s3d_scene_view* view_samp = NULL;
int has_sampled = 0;
- int has_receiver = 0;
+ int cpt_receiver = 0;
res_T res = RES_OK;
ASSERT(scn && out_view_rt && out_view_samp);
@@ -297,8 +297,21 @@ scene_create_s3d_views
unsigned id;
htable_instance_iterator_next(&it);
- if(inst->receiver_mask) {
- has_receiver = 1;
+ if (inst->receiver_mask & SSOL_FRONT) {
+ ASSERT(SSOL_FRONT == 1 || SSOL_FRONT == 2);
+ if (cpt_receiver == INT_MAX) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ inst->mc_result_idx[side_idx(SSOL_FRONT)] = cpt_receiver++;
+ }
+ if (inst->receiver_mask & SSOL_BACK) {
+ ASSERT(SSOL_BACK == 1 || SSOL_BACK == 2);
+ if (cpt_receiver == INT_MAX) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ inst->mc_result_idx[side_idx(SSOL_BACK)] = cpt_receiver++;
}
if(!inst->sample) continue;
@@ -327,7 +340,7 @@ scene_create_s3d_views
goto error;
}
- if(!has_receiver) {
+ if(!cpt_receiver) {
log_warning(scn->dev, "No receiver is defined.\n");
}
@@ -337,6 +350,7 @@ scene_create_s3d_views
if(res != RES_OK) goto error;
exit:
+ scn->nb_receivers = cpt_receiver;
*out_view_rt = view_rt;
*out_view_samp = view_samp;
return res;
diff --git a/src/ssol_scene_c.h b/src/ssol_scene_c.h
@@ -36,7 +36,7 @@ struct ssol_scene;
struct ssol_sun;
struct ssol_scene {
- /* Map the instantiated RT/Samp S3D shape id to its SSOL intrance */
+ /* Map the instantiated RT/Samp S3D shape id to its SSOL instance */
struct htable_instance instances_rt;
struct htable_instance instances_samp;
@@ -46,6 +46,8 @@ struct ssol_scene {
struct ssol_sun* sun; /* Sun of the scene */
struct ssol_atmosphere* atmosphere; /* Atmosphere of the scene */
+ int nb_receivers;
+
struct ssol_device* dev;
ref_T ref;
};
diff --git a/src/ssol_solver.c b/src/ssol_solver.c
@@ -301,6 +301,15 @@ check_scene(const struct ssol_scene* scene, const char* caller)
return RES_OK;
}
+static FINLINE int
+get_rcv_idx(const struct point* pt)
+{
+ ASSERT(pt);
+ if (pt->inst->receiver_mask & pt->side)
+ return pt->inst->mc_result_idx[side_idx(pt->side)];
+ return MC_RESULT_IDX_NONE;
+}
+
/*******************************************************************************
* Exported functions
******************************************************************************/
@@ -321,10 +330,12 @@ ssol_solve
struct ssp_rng_proxy* rng_proxy = NULL;
struct mc_data* mc_shadows = NULL;
struct mc_data* mc_missings = NULL;
+ struct mc_data* mc_global_receivers = NULL;
float sampled_area = 0;
int nthreads = 0;
int nrealisations = 0;
- int i = 0;
+ int nb_receivers = 0;
+ int i = 0, j =0;
ATOMIC res = RES_OK;
ASSERT(nrealisations <= INT_MAX);
@@ -341,6 +352,7 @@ ssol_solve
goto error;
}
nrealisations = (int)realisations_count;
+ nthreads = (int) scn->dev->nthreads;
res = check_scene(scn, FUNC_NAME);
if(res != RES_OK) goto error;
@@ -352,28 +364,35 @@ ssol_solve
if(res != RES_OK) goto error;
S3D(scene_view_compute_area(view_samp, &sampled_area));
+ /* allocate room for per-receiver global MC data */
+ nb_receivers = scn->nb_receivers;
+ res = darray_receiver_resize(&estimator->global_receiver, nb_receivers);
+ if (res != RES_OK) goto error;
+
/* Create a RNG proxy from the submitted RNG state */
res = ssp_rng_proxy_create_from_rng
- (scn->dev->allocator, rng_state, scn->dev->nthreads, &rng_proxy);
+ (scn->dev->allocator, rng_state, nthreads, &rng_proxy);
if(res != RES_OK) goto error;
/* Create per thread data structures */
- #define CREATE(Data) { \
+ #define CREATEN(Data,N) { \
ASSERT(!(Data)); \
- if(!sa_add((Data), scn->dev->nthreads)) { \
+ if(!sa_add((Data), nthreads)) { \
res = RES_BAD_ARG; \
goto error; \
} \
- memset((Data), 0, sizeof((Data)[0])*scn->dev->nthreads); \
+ memset((Data), 0, sizeof((Data)[0])*N*nthreads); \
} (void)0
+ #define CREATE(Data) CREATEN(Data,1)
CREATE(rngs);
CREATE(bsdfs);
CREATE(mc_shadows);
CREATE(mc_missings);
+ CREATEN(mc_global_receivers, nb_receivers);
#undef CREATE
+ #undef CREATEN
/* Setup per thread data structures */
- nthreads = (int)scn->dev->nthreads;
FOR_EACH(i, 0, nthreads) {
res = ssf_bsdf_create(scn->dev->allocator, bsdfs+i);
if(res != RES_OK) goto error;
@@ -389,6 +408,7 @@ ssol_solve
struct ssf_bsdf* bsdf;
struct mc_data* shadow;
struct mc_data* missing;
+ struct mc_data* global_receiver;
float org[3], dir[3], range[2] = { 0, FLT_MAX };
const int ithread = omp_get_thread_num();
size_t depth = 0;
@@ -402,6 +422,7 @@ ssol_solve
bsdf = bsdfs[ithread];
shadow = mc_shadows + ithread;
missing = mc_missings + ithread;
+ global_receiver = mc_global_receivers + ithread * nb_receivers;
/* Find a new starting point of the radiative random walk */
is_lit = point_init(&pt, scn, sampled_area, view_samp, view_rt,
@@ -425,10 +446,21 @@ ssol_solve
if(point_is_receiver(&pt)) {
const res_T res_local = point_dump(&pt, (size_t)i, depth, output);
+ int rcv_idx = MC_RESULT_IDX_NONE;
+ struct mc_data* global_recv = NULL;
if(res_local != RES_OK) {
ATOMIC_SET(&res, res_local);
break;
}
+ rcv_idx = get_rcv_idx(&pt);
+ if (rcv_idx == MC_RESULT_IDX_NONE) {
+ ATOMIC_SET(&res, res_local);
+ break;
+ }
+ ASSERT(rcv_idx < nb_receivers);
+ global_recv = &global_receiver[rcv_idx];
+ global_recv->weight += pt.weight;
+ global_recv->sqr_weight += pt.weight*pt.weight;
hit_a_receiver = 1;
}
@@ -484,11 +516,19 @@ ssol_solve
}
/* Merge per thread estimations */
+ ASSERT(estimator->global_receiver.size <= nb_receivers);
FOR_EACH(i, 0, nthreads) {
estimator->shadow.weight += mc_shadows[i].weight;
estimator->shadow.sqr_weight += mc_shadows[i].sqr_weight;
estimator->missing.weight += mc_missings[i].weight;
estimator->missing.sqr_weight += mc_missings[i].sqr_weight;
+ FOR_EACH(j, 0, nb_receivers) {
+ const int idx = i * nb_receivers + j;
+ estimator->global_receiver.data[j].weight
+ += mc_global_receivers[idx].weight;
+ estimator->global_receiver.data[j].sqr_weight
+ += mc_global_receivers[idx].sqr_weight;
+ }
}
estimator->realisation_count += realisations_count;
@@ -508,6 +548,7 @@ exit:
}
sa_release(mc_shadows);
sa_release(mc_missings);
+ sa_release(mc_global_receivers);
return (res_T)res;
error:
goto exit;
diff --git a/src/test_ssol_estimator.c b/src/test_ssol_estimator.c
@@ -38,6 +38,10 @@ main(int argc, char** argv)
struct ssol_device* dev;
struct ssp_rng* rng;
struct ssol_estimator* estimator;
+ struct ssol_instance* inst;
+ struct ssol_material* v_mtl;
+ struct ssol_shape* shape;
+ struct ssol_object* object;
struct ssol_estimator_status status;
size_t count;
@@ -71,6 +75,21 @@ main(int argc, char** argv)
CHECK(GET_STATUS(estimator, SSOL_STATUS_MISSING, &status), RES_OK);
#undef GET_STATUS
+ #define GET_RCV_STATUS ssol_estimator_get_receiver_status
+ CHECK(ssol_material_create_virtual(dev, &v_mtl), RES_OK);
+ CHECK(ssol_object_create(dev, &object), RES_OK);
+ CHECK(ssol_shape_create_punched_surface(dev, &shape), RES_OK);
+ CHECK(ssol_object_add_shaded_shape(object, shape, v_mtl, v_mtl), RES_OK);
+ CHECK(ssol_object_instantiate(object, &inst), RES_OK);
+ CHECK(ssol_instance_set_receiver(inst, SSOL_FRONT), RES_OK);
+ CHECK(GET_RCV_STATUS(NULL, inst, SSOL_BACK, &status), RES_BAD_ARG);
+ CHECK(GET_RCV_STATUS(estimator, NULL, SSOL_BACK, &status), RES_BAD_ARG);
+ CHECK(GET_RCV_STATUS(estimator, inst, 0, &status), RES_BAD_ARG);
+ CHECK(GET_RCV_STATUS(estimator, inst, SSOL_BACK, NULL), RES_BAD_ARG);
+ /* we cannot check that a status is available for the front face
+ solve has been succesfully called */
+ #undef GET_RCV_STATUS
+
CHECK(ssol_estimator_get_count(NULL, &count), RES_BAD_ARG);
CHECK(ssol_estimator_get_count(estimator, NULL), RES_BAD_ARG);
CHECK(ssol_estimator_get_count(estimator, &count), RES_OK);
@@ -79,6 +98,10 @@ main(int argc, char** argv)
CHECK(ssol_estimator_get_failed_count(estimator, NULL), RES_BAD_ARG);
CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK);
+ CHECK(ssol_material_ref_put(v_mtl), RES_OK);
+ CHECK(ssol_object_ref_put(object), RES_OK);
+ CHECK(ssol_shape_ref_put(shape), RES_OK);
+ CHECK(ssol_instance_ref_put(inst), RES_OK);
CHECK(ssol_estimator_ref_put(estimator), RES_OK);
CHECK(ssol_device_ref_put(dev), RES_OK);
diff --git a/src/test_ssol_solver2.c b/src/test_ssol_solver2.c
@@ -191,13 +191,30 @@ main(int argc, char** argv)
CHECK(eq_eps(m, 4 * DNI_cos, 4 * DNI_cos * 1e-4), 1);
#define SQR(x) ((x)*(x))
CHECK(eq_eps(std, 0, 1e-4), 1);
- CHECK(ssol_estimator_get_status(estimator, SSOL_STATUS_SHADOW, &status), RES_OK);
+#define GET_STATUS ssol_estimator_get_status
+ CHECK(GET_STATUS(estimator, SSOL_STATUS_SHADOW, &status), RES_OK);
logger_print(&logger, LOG_OUTPUT, "Shadows = %g +/- %g", status.E, status.SE);
CHECK(eq_eps(status.E, 0, 1e-4), 1);
- CHECK(ssol_estimator_get_status(estimator, SSOL_STATUS_MISSING, &status), RES_OK);
+ CHECK(GET_STATUS(estimator, SSOL_STATUS_MISSING, &status), RES_OK);
logger_print(&logger, LOG_OUTPUT, "Missing = %g +/- %g", status.E, status.SE);
CHECK(eq_eps(status.E, 0, 1e-4), 1);
+#define GET_RCV_STATUS ssol_estimator_get_receiver_status
+ CHECK(GET_RCV_STATUS(estimator, heliostat1, SSOL_BACK, &status), RES_BAD_ARG);
+ CHECK(GET_RCV_STATUS(estimator, heliostat1, SSOL_FRONT, &status), RES_OK);
+ logger_print(&logger, LOG_OUTPUT, "P(heliostat1) = %g +/- %g", status.E, status.SE);
+ CHECK(GET_RCV_STATUS(estimator, heliostat2, SSOL_BACK, &status), RES_BAD_ARG);
+ CHECK(GET_RCV_STATUS(estimator, heliostat2, SSOL_FRONT, &status), RES_OK);
+ logger_print(&logger, LOG_OUTPUT, "P(heliostat2) = %g +/- %g", status.E, status.SE);
+ CHECK(GET_RCV_STATUS(estimator, secondary, SSOL_BACK, &status), RES_BAD_ARG);
+ CHECK(GET_RCV_STATUS(estimator, secondary, SSOL_FRONT, &status), RES_OK);
+ logger_print(&logger, LOG_OUTPUT, "P(secondary) = %g +/- %g", status.E, status.SE);
+ CHECK(GET_RCV_STATUS(estimator, target, SSOL_BACK, &status), RES_BAD_ARG);
+ CHECK(GET_RCV_STATUS(estimator, target, SSOL_FRONT, &status), RES_OK);
+ logger_print(&logger, LOG_OUTPUT, "P(target) = %g +/- %g", status.E, status.SE);
+ CHECK(eq_eps(status.E, m, 1e-8), 1);
+ CHECK(eq_eps(status.SE, std, 1e-4), 1);
+
/* Free data */
CHECK(ssol_instance_ref_put(heliostat1), RES_OK);
CHECK(ssol_instance_ref_put(heliostat2), RES_OK);
