commit de44f0c3189f33c8115036c6f2d9fa51999e372c
parent c116b67b0b45168552d0f6a9b9af4ca335a2510e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 25 Sep 2015 16:32:36 +0200
Add and test the sgf_estimator_get_status_infinity function
Return the estimated radiative flux that does not reach a geometry.
Diffstat:
4 files changed, 105 insertions(+), 16 deletions(-)
diff --git a/src/sgf.h b/src/sgf.h
@@ -131,6 +131,8 @@ SGF_API res_T
sgf_estimator_ref_put
(struct sgf_estimator* estimator);
+/* Return the estimated Gebhart Factor between the estimator primitive and
+ * `primitive_id'. */
SGF_API res_T
sgf_estimator_get_status
(struct sgf_estimator* estimator,
@@ -138,6 +140,13 @@ sgf_estimator_get_status
const size_t spectral_band_id,
struct sgf_status* status);
+/* Return the estimated radiative flux that hit nothing */
+SGF_API res_T
+sgf_estimator_get_status_infinity
+ (struct sgf_estimator* estimator,
+ const size_t spectral_band,
+ struct sgf_status* status);
+
END_DECLS
#endif /* SGF_H */
diff --git a/src/sgf_estimator.c b/src/sgf_estimator.c
@@ -49,6 +49,8 @@ struct accum {
/* Estimator of the Gebhart Factor of a given face to all the other ones */
struct sgf_estimator {
struct darray_accum buf; /* Integrated radiative flux of each primitive */
+ struct darray_accum buf_infinity; /* Radiative flux going to the infinity */
+
size_t nsteps; /* # realisations */
size_t nprims; /* # primitives */
size_t nbands; /* # spectral bands */
@@ -75,10 +77,27 @@ accum_weight(struct accum* accums, const size_t iface, const double weight)
accums[iface].sqr_radiative_flux += weight * weight;
}
+static FINLINE void
+setup_status
+ (const struct accum* acc,
+ const size_t nsteps,
+ struct sgf_status* status)
+{
+ ASSERT(acc && status && nsteps);
+
+ status->E = acc->radiative_flux / (double)nsteps; /* Expected value */
+ status->V = /* Variance */
+ acc->sqr_radiative_flux / (double)nsteps
+ - (acc->radiative_flux * acc->radiative_flux) / (double)(nsteps * nsteps);
+ status->SE = sqrt(status->V / (double)nsteps); /* Standard error */
+ status->nsteps = nsteps; /* # realisations */
+}
+
/* Return RES_BAD_OP on numerical issue. i.e. the radiative path is skipped */
static res_T
gebhart_radiative_path
(struct accum* accums,
+ struct accum* accum_infinity,
struct ssp_rng* rng,
const size_t primitive_id,
struct sgf_scene_desc* desc)
@@ -88,10 +107,10 @@ gebhart_radiative_path
struct s3d_attrib attrib;
const double trans_min = 1.e-8; /* Minimum transmissivity threshold */
- double outside_radiative_flux = 0.0;
double proba_reflec_spec;
double transmissivity, emissivity, reflectivity, specularity;
#ifndef NDEBUG
+ double infinite_radiative_flux = 0.0;
double sum_radiative_flux = 0.f;
#endif
@@ -103,6 +122,8 @@ gebhart_radiative_path
float range[2]; /* Traced ray range */
float u, v;
+ ASSERT(accums && accum_infinity && rng && desc);
+
/* Discard faces with no emissivity */
emissivity = desc->get_material_property(desc->material,
SGF_MATERIAL_EMISSIVITY, primitive_id, 0);
@@ -127,7 +148,6 @@ gebhart_radiative_path
ssp_ran_hemisphere_cos(rng, normal, dir);
transmissivity = 1.0;
- outside_radiative_flux = 0.0;
for(;;) { /* Here we go */
prim_from = prim;
@@ -140,7 +160,11 @@ gebhart_radiative_path
}
if(S3D_HIT_NONE(&hit)) { /* The ray is outside the volume */
- outside_radiative_flux += transmissivity;
+ accum_infinity->radiative_flux += transmissivity;
+ accum_infinity->sqr_radiative_flux += transmissivity * transmissivity;
+#ifndef NDEBUG
+ infinite_radiative_flux = transmissivity;
+#endif
break;
}
@@ -194,7 +218,7 @@ gebhart_radiative_path
}
#if !defined(NDEBUG)
/* Check the energy conservation property */
- ASSERT(eq_eps(sum_radiative_flux + outside_radiative_flux, 1.0, 1.e6));
+ ASSERT(eq_eps(sum_radiative_flux + infinite_radiative_flux, 1.0, 1.e6));
#endif
return RES_OK;
}
@@ -219,6 +243,7 @@ estimator_create(struct sgf_device* dev, struct sgf_estimator** out_estimator)
SGF(device_ref_get(dev));
estimator->dev = dev;
darray_accum_init(dev->allocator, &estimator->buf);
+ darray_accum_init(dev->allocator, &estimator->buf_infinity);
exit:
*out_estimator = estimator;
@@ -240,6 +265,7 @@ estimator_release(ref_T* ref)
ASSERT(ref);
estimator = CONTAINER_OF(ref, struct sgf_estimator, ref);
darray_accum_release(&estimator->buf);
+ darray_accum_release(&estimator->buf_infinity);
dev = estimator->dev;
MEM_RM(dev->allocator, estimator);
SGF(device_ref_put(dev));
@@ -295,24 +321,40 @@ sgf_integrate
goto error;
}
- /* Allocate internal memory space */
+ /* Allocate and init the Gebhart Factor result buffer */
res = darray_accum_resize(&estimator->buf, nprims*desc->spectral_bands_count);
if(res != RES_OK) {
if(dev->verbose) {
logger_print(dev->logger, LOG_ERROR,
-"Couldn't allocate the Gebhart Factor integration data structure.\n");
+ "Couldn't allocate the Gebhart Factor result buffer.\n");
}
goto error;
}
memset(darray_accum_data_get(&estimator->buf), 0,
darray_accum_size_get(&estimator->buf)*sizeof(struct accum));
+ /* Allocate and init the infinite radiative flux buffer */
+ res = darray_accum_resize(&estimator->buf_infinity, desc->spectral_bands_count);
+ if(res != RES_OK) {
+ if(dev->verbose) {
+ logger_print(dev->logger, LOG_ERROR,
+ "Couldn't allocate the buffer of infinite radiative flux.\n");
+ }
+ goto error;
+ }
+ memset(darray_accum_data_get(&estimator->buf_infinity), 0,
+ darray_accum_size_get(&estimator->buf_infinity)*sizeof(struct accum));
+
/* Invoke the Gebhart factor integration. */
FOR_EACH(ispectral_band, 0, desc->spectral_bands_count) {
struct accum* accums =
darray_accum_data_get(&estimator->buf) + ispectral_band * nprims;
+ struct accum* accum_infinity =
+ darray_accum_data_get(&estimator->buf_infinity) + ispectral_band;
+
FOR_EACH(istep, 0, steps_count) {
- res = gebhart_radiative_path(accums, rng, primitive_id, desc);
+ res = gebhart_radiative_path
+ (accums, accum_infinity, rng, primitive_id, desc);
if(res != RES_OK) goto error;
}
}
@@ -352,7 +394,6 @@ sgf_estimator_get_status
struct sgf_status* status)
{
struct accum* acc;
- size_t nsteps;
size_t iacc;
if(!estimator || !status)
@@ -376,14 +417,33 @@ sgf_estimator_get_status
iacc = ispectral_band * estimator->nprims + iprimitive;
ASSERT(iacc < darray_accum_size_get(&estimator->buf));
acc = darray_accum_data_get(&estimator->buf) + iacc;
- nsteps = estimator->nsteps;
+ setup_status(acc, estimator->nsteps, status);
- status->E = acc->radiative_flux / (double)nsteps; /* Expected value */
- status->V = /* Variance */
- acc->sqr_radiative_flux / (double)nsteps
- - (acc->radiative_flux * acc->radiative_flux) / (double)(nsteps * nsteps);
- status->SE = sqrt(status->V / (double)nsteps); /* Standard error */
- status->nsteps = nsteps; /* # realisations */
+ return RES_OK;
+}
+
+res_T
+sgf_estimator_get_status_infinity
+ (struct sgf_estimator* estimator,
+ const size_t ispectral_band,
+ struct sgf_status* status)
+{
+ struct accum* acc;
+
+ if(!estimator || !status)
+ return RES_BAD_ARG;
+
+ if(ispectral_band >= estimator->nbands) {
+ if(estimator->dev->verbose) {
+ logger_print(estimator->dev->logger, LOG_ERROR,
+ "Out of bound spectral band index `%lu'\n", ispectral_band);
+ }
+ return RES_BAD_ARG;
+ }
+
+ ASSERT(ispectral_band < darray_accum_size_get(&estimator->buf_infinity));
+ acc = darray_accum_data_get(&estimator->buf_infinity) + ispectral_band;
+ setup_status(acc, estimator->nsteps, status);
return RES_OK;
}
diff --git a/src/test_sgf_cube.c b/src/test_sgf_cube.c
@@ -232,7 +232,17 @@ main(int argc, char** argv)
CHECK(sgf_estimator_get_status(estimator, SIZE_MAX, 0, status), RES_BAD_ARG);
CHECK(sgf_estimator_get_status(NULL, 0, 0, status), RES_BAD_ARG);
CHECK(sgf_estimator_get_status(estimator, 0, 0, status), RES_OK);
+ CHECK(status->nsteps, NSTEPS);
+ CHECK(sgf_estimator_get_status_infinity(NULL, SIZE_MAX, NULL), RES_BAD_ARG);
+ CHECK(sgf_estimator_get_status_infinity(estimator, SIZE_MAX, NULL), RES_BAD_ARG);
+ CHECK(sgf_estimator_get_status_infinity(NULL, 0, NULL), RES_BAD_ARG);
+ CHECK(sgf_estimator_get_status_infinity(estimator, 0, NULL), RES_BAD_ARG);
+ CHECK(sgf_estimator_get_status_infinity(NULL, SIZE_MAX, NULL), RES_BAD_ARG);
+ CHECK(sgf_estimator_get_status_infinity(estimator, SIZE_MAX, status), RES_BAD_ARG);
+ CHECK(sgf_estimator_get_status_infinity(NULL, 0, status), RES_BAD_ARG);
+ CHECK(sgf_estimator_get_status_infinity(estimator, 0, status), RES_OK);
+ CHECK(eq_eps(status->E, 0, status->SE), 1);
CHECK(status->nsteps, NSTEPS);
CHECK(sgf_estimator_ref_get(NULL), RES_BAD_ARG);
@@ -246,15 +256,20 @@ main(int argc, char** argv)
for(iprim = 0; iprim < nprims; ++iprim) {
size_t iprim2;
struct sgf_status* row = status + iprim * nprims;
+ struct sgf_status inf;
const int ithread = omp_get_thread_num();
struct sgf_estimator* est;
CHECK(sgf_integrate(sgf, rngs[ithread], iprim, &scn_desc, NSTEPS, &est), RES_OK);
+ CHECK(sgf_estimator_get_status_infinity(est, 0, &inf), RES_OK);
+ CHECK(eq_eps(inf.E, 0, 3 * inf.SE), 1);
+
FOR_EACH(iprim2, 0, nprims) {
CHECK(sgf_estimator_get_status(est, iprim2, 0, row + iprim2), RES_OK);
CHECK(row[iprim2].nsteps, NSTEPS);
}
+
CHECK(sgf_estimator_ref_put(est), RES_OK);
}
diff --git a/src/test_sgf_tetrahedron.c b/src/test_sgf_tetrahedron.c
@@ -121,11 +121,16 @@ main(int argc, char** argv)
#pragma omp parallel for
for(iprim = 0; iprim < nprims; ++iprim) {
struct sgf_status* row = status + iprim*nprims;
+ struct sgf_status inf;
struct sgf_estimator* estimator = NULL;
size_t iprim2;
const int ithread = omp_get_thread_num();
CHECK(sgf_integrate(sgf, rngs[ithread], iprim, &desc, NSTEPS, &estimator), RES_OK);
+
+ CHECK(sgf_estimator_get_status_infinity(estimator, 0, &inf), RES_OK);
+ CHECK(eq_eps(inf.E, 0, 3*inf.SE), 1);
+
FOR_EACH(iprim2, 0, nprims) {
CHECK(sgf_estimator_get_status(estimator, iprim2, 0, row + iprim2), RES_OK);
CHECK(row[iprim2].nsteps, NSTEPS);
@@ -145,7 +150,7 @@ main(int argc, char** argv)
}
printf("\n");
/* Ensure the energy conservation property */
- CHECK(eq_eps(sum, 1.0, 1.e-6), 1);
+ CHECK(eq_eps(sum, 1.0, 1.e-3), 1);
}
printf("\n");
