commit 55800a0a35915ba325d91bc8f87b97af61b154a4
parent 6901f6b770584b2175eae2e8e446c3f335693176
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Sat, 16 Sep 2017 15:11:40 +0200
Merge branch 'release_0.4'
Diffstat:
8 files changed, 182 insertions(+), 69 deletions(-)
diff --git a/README.md b/README.md
@@ -26,6 +26,15 @@ variable the install directories of its dependencies.
## Release notes
+### Version 0.4
+
+- Add the `SSOL_PATH_ERROR` type used for the paths that travel unforeseen
+ mediums.
+- Fix the cosine factor estimation that did not take into account the
+ shadowed realisations.
+- Ensure the energy conservation property for dielectric materials. Previously,
+ some energy was lost even for dielectric materials with no absorption.
+
### Version 0.3
- Full rewrite of the estimated values. The global results report the cosine
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -28,7 +28,7 @@ find_package(RSys 0.4 REQUIRED)
find_package(Star3D 0.4.1 REQUIRED)
find_package(Star3DUT 0.2 REQUIRED)
find_package(StarCPR 0.1 REQUIRED)
-find_package(StarSF 0.2 REQUIRED)
+find_package(StarSF 0.1 REQUIRED)
find_package(StarSP 0.4 REQUIRED)
find_package(OpenMP 1.2 REQUIRED)
@@ -50,7 +50,7 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D Star3DUT StarCPR StarSF Sta
# Configure and define targets
################################################################################
set(VERSION_MAJOR 0)
-set(VERSION_MINOR 3)
+set(VERSION_MINOR 4)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
@@ -130,8 +130,20 @@ rcmake_setup_devel(solstice-solver SolSolver ${VERSION} solstice/ssol_version.h)
# Add tests
################################################################################
if(NOT NO_TEST)
+ set(SSOL_TEST_FILES_INC
+ test_ssol_circ2D_geometry.h
+ test_ssol_geometries.h
+ test_ssol_rect_geometry.h
+ test_ssol_utils.h
+ test_ssol_cube_geometry.h
+ test_ssol_materials.h
+ test_ssol_rect2D_geometry.h
+ )
+
+ rcmake_prepend_path(SSOL_TEST_FILES_INC ${SSOL_SOURCE_DIR})
+
function(build_test _name)
- add_executable(${_name} ${SSOL_SOURCE_DIR}/${_name}.c)
+ add_executable(${_name} ${SSOL_SOURCE_DIR}/${_name}.c ${SSOL_TEST_FILES_INC})
target_link_libraries(${_name}
solstice-solver RSys Star3D StarSP)
endfunction()
@@ -174,7 +186,6 @@ if(NOT NO_TEST)
build_test(test_ssol_draw)
register_test(test_ssol_draw_draft test_ssol_draw draft)
register_test(test_ssol_draw_pt test_ssol_draw pt)
-
endif()
################################################################################
diff --git a/src/ssol.h b/src/ssol.h
@@ -68,8 +68,9 @@ enum ssol_side_flag {
enum ssol_path_type {
SSOL_PATH_MISSING, /* The path misses the receivers */
- SSOL_PATH_SHADOW, /* The path is occluded before the sampled geometry */
- SSOL_PATH_SUCCESS /* The path contributes to at least one receiver */
+ SSOL_PATH_SHADOW, /* The path is occluded before the sampled geometry */
+ SSOL_PATH_SUCCESS, /* The path contributes to at least one receiver */
+ SSOL_PATH_ERROR /* The path was canceled due to a path-related error */
};
enum ssol_material_type {
diff --git a/src/ssol_estimator_c.h b/src/ssol_estimator_c.h
@@ -101,6 +101,14 @@ mc_data_get(struct mc_data* data, double* weight, double* sqr_weight)
*sqr_weight = data->sqr_weight__;
}
+static INLINE void
+mc_data_cancel(struct mc_data* data, size_t irealisation)
+{
+ ASSERT(data);
+ if(data->irealisation!=irealisation) return;
+ data->tmp = 0;
+}
+
/*******************************************************************************
* One sided per shape MC data
******************************************************************************/
diff --git a/src/ssol_material.c b/src/ssol_material.c
@@ -85,9 +85,7 @@ setup_dielectric_bsdf
const struct ssol_medium* medium,
struct ssf_bsdf* bsdf)
{
- struct ssf_bxdf* brdf = NULL;
- struct ssf_bxdf* btdf = NULL;
- struct ssf_fresnel* fresnel = NULL;
+ struct ssf_bxdf* bxdf = NULL;
double eta_i, eta_t;
res_T res = RES_OK;
ASSERT(mtl && fragment && mtl->type == SSOL_MATERIAL_DIELECTRIC);
@@ -103,25 +101,19 @@ setup_dielectric_bsdf
eta_i = ssol_data_get_value(&mtl->out_medium.refractive_index, wavelength);
eta_t = ssol_data_get_value(&mtl->in_medium.refractive_index, wavelength);
- #define CALL(Func) { res = Func; if(res != RES_OK) goto error; } (void)0
- /* Setup the reflective part */
- CALL(ssf_fresnel_create
- (mtl->dev->allocator, &ssf_fresnel_dielectric_dielectric, &fresnel));
- CALL(ssf_fresnel_dielectric_dielectric_setup(fresnel, eta_i, eta_t));
- CALL(ssf_bxdf_create(mtl->dev->allocator, &ssf_specular_reflection, &brdf));
- CALL(ssf_specular_reflection_setup(brdf, fresnel));
- /* Setup the transmissive part */
- CALL(ssf_bxdf_create(mtl->dev->allocator, &ssf_specular_transmission, &btdf));
- CALL(ssf_specular_transmission_setup(btdf, eta_i, eta_t));
- /* Setup the scattering function */
- CALL(ssf_bsdf_add(bsdf, brdf, 0.5));
- CALL(ssf_bsdf_add(bsdf, btdf, 0.5));
- #undef CALL
+ res = ssf_bxdf_create
+ (mtl->dev->allocator, &ssf_specular_dielectric_dielectric_interface, &bxdf);
+ if(res != RES_OK) goto error;
+
+ res = ssf_specular_dielectric_dielectric_interface_setup
+ (bxdf, eta_i, eta_t);
+ if(res != RES_OK) goto error;
+
+ res = ssf_bsdf_add(bsdf, bxdf, 1.0);
+ if(res != RES_OK) goto error;
exit:
- if(brdf) SSF(bxdf_ref_put(brdf));
- if(btdf) SSF(bxdf_ref_put(btdf));
- if(fresnel) SSF(fresnel_ref_put(fresnel));
+ if(bxdf) SSF(bxdf_ref_put(bxdf));
return res;
error:
goto exit;
diff --git a/src/ssol_solver.c b/src/ssol_solver.c
@@ -173,7 +173,8 @@ struct point {
const struct ssol_material* material;
/* tmp quantities to compute weights */
double kabs_at_pt;
- double partial_atm, partial_recv, partial_other;
+ /* for conservation of energy check */
+ double energy_loss;
/* MC weights */
/* Set once */
double initial_flux; /* the initial flux*/
@@ -204,7 +205,8 @@ struct point {
{0, 0}, /* UV */ \
0, /* Wavelength */ \
NULL, /* Material */ \
- 0, 0, 0, 0, /* tmp values */ \
+ 0, /* tmp values */ \
+ 0, /* Energy loss */ \
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, /* MC weights */ \
SSOL_FRONT /* Side */ \
}
@@ -306,7 +308,7 @@ point_init
pt->cos_factor = surface_sun_cos;
}
- pt->partial_atm = pt->partial_recv = pt->partial_other = 0;
+ pt->energy_loss = w0;
pt->initial_flux = w0;
pt->prev_outgoing_flux = w0;
pt->prev_outgoing_if_no_atm_loss = w0;
@@ -482,7 +484,7 @@ point_shade
pt->kabs_at_pt = (1 - propagated);
pt->outgoing_flux = pt->incoming_flux * propagated;
pt->outgoing_if_no_atm_loss = pt->incoming_if_no_atm_loss * propagated;
- pt->outgoing_if_no_field_loss = point_is_receiver(pt)
+ pt->outgoing_if_no_field_loss = point_is_receiver(pt)
? pt->incoming_if_no_field_loss*propagated : pt->incoming_if_no_field_loss;
if(type & SSF_TRANSMISSION) {
@@ -665,12 +667,6 @@ error:
goto exit;
}
-/*
- * FIXME Are the following accumulations OK when the radiative path bounces
- * several times on the same receiver ? It seems weird to add the overall
- * absorptivity and reflectivity losses to the corresponding per receiver
- * accumulators since they already registered some losses.
- */
static res_T
update_mc
(struct point* pt,
@@ -682,7 +678,11 @@ update_mc
res_T res = RES_OK;
ASSERT(pt && thread_ctx && point_is_receiver(pt));
- pt->partial_recv += pt->incoming_flux - pt->outgoing_flux;
+ #define ACCUM_WEIGHT(Name, W)\
+ mc_data_add_weight(&thread_ctx->Name, irealisation, W)
+ ACCUM_WEIGHT(absorbed_by_receivers, pt->incoming_flux - pt->outgoing_flux);
+ pt->energy_loss -= (pt->incoming_flux - pt->outgoing_flux);
+ #undef ACCUM_WEIGHT
/* Per receiver MC accumulation */
res = get_mc_receiver_1side(&thread_ctx->mc_rcvs, pt->inst, pt->side, &mc_rcv1);
@@ -774,7 +774,7 @@ trace_radiative_path
/* Find a new starting point of the radiative random walk */
res = point_init(&pt, sampled_area_proxy, scn, &thread_ctx->mc_samps,
- view_samp, view_rt, ran_sun_dir, ran_sun_wl, thread_ctx->rng,
+ view_samp, view_rt, ran_sun_dir, ran_sun_wl, thread_ctx->rng,
&in_medium, &is_lit);
if(res != RES_OK) goto error;
@@ -799,6 +799,7 @@ trace_radiative_path
if(!is_lit) { /* The starting point is not lit */
ACCUM_WEIGHT(pt.mc_samp->shadowed, pt.initial_flux);
ACCUM_WEIGHT(thread_ctx->shadowed, pt.initial_flux);
+ pt.energy_loss -= pt.initial_flux;
if(tracker) path.type = SSOL_PATH_SHADOW;
} else {
/* Setup the ray as if it starts from the current point position in order
@@ -846,7 +847,9 @@ trace_radiative_path
res = update_mc(&pt, irealisation, thread_ctx);
if(res != RES_OK) goto error;
} else {
- pt.partial_other += pt.incoming_flux * pt.kabs_at_pt;
+ ACCUM_WEIGHT(thread_ctx->other_absorbed,
+ pt.incoming_flux * pt.kabs_at_pt);
+ pt.energy_loss -= (pt.incoming_flux * pt.kabs_at_pt);
}
/* Stop the radiative random walk if no more flux */
@@ -888,7 +891,11 @@ trace_radiative_path
if (res != RES_OK) goto error;
}
last_segment = 1; /* Path reached its last segment */
- ASSERT(in_atm);
+ if(!in_atm) {
+ log_error(scn->dev, "Inconsistent medium description.\n");
+ res = RES_BAD_OP;
+ goto error;
+ }
}
/* Don't change prev_outgoing weigths nor record segment absorption until
@@ -897,10 +904,14 @@ trace_radiative_path
* a non-virtual material is hit or no further hit can be found. */
if(last_segment || !hit_virtual) {
if(in_atm) {
- pt.partial_atm += pt.prev_outgoing_flux - pt.incoming_flux;
+ ACCUM_WEIGHT(thread_ctx->absorbed_by_atmosphere,
+ pt.prev_outgoing_flux - pt.incoming_flux);
} else {
- pt.partial_other += pt.prev_outgoing_flux - pt.incoming_flux;
+ ACCUM_WEIGHT(thread_ctx->other_absorbed,
+ pt.prev_outgoing_flux - pt.incoming_flux);
}
+ pt.energy_loss -= (pt.prev_outgoing_flux - pt.incoming_flux);
+
if(last_segment) {
break;
}
@@ -910,7 +921,8 @@ trace_radiative_path
}
depth += !hit_virtual;
- ASSERT(depth < 100); /* FIXME: create a true cancel path for this MC sample */
+ /* FIXME: create a true cancel path for this MC sample */
+ ASSERT(depth < 100);
/* Update the point */
point_update_from_hit(&pt, scn, org, dir, &hit, &ray_data);
@@ -919,42 +931,122 @@ trace_radiative_path
res = path_add_vertex(&path, pt.pos, pt.outgoing_flux);
if (res != RES_OK) goto error;
}
-
+
ssol_medium_copy(&in_medium, &out_medium);
}
- /* Check conservation of energy */
- ASSERT((double)depth * pt.initial_flux * DBL_EPSILON >=
- fabs(pt.initial_flux -
- (pt.outgoing_flux + pt.partial_recv + pt.partial_atm + pt.partial_other)));
-
- /* Now that the sample ends successfully,
- * record MC weights and register the remaining power as missing */
- ACCUM_WEIGHT(pt.mc_samp->cos_factor, pt.cos_factor);
- ACCUM_WEIGHT(thread_ctx->cos_factor, pt.cos_factor);
+
+ /* Register the remaining flux as missing */
ACCUM_WEIGHT(thread_ctx->missing, pt.outgoing_flux);
- ACCUM_WEIGHT(thread_ctx->absorbed_by_receivers, pt.partial_recv);
- ACCUM_WEIGHT(thread_ctx->absorbed_by_atmosphere, pt.partial_atm);
- ACCUM_WEIGHT(thread_ctx->other_absorbed, pt.partial_other);
- #undef ACCUM_WEIGHT
+ pt.energy_loss -= pt.outgoing_flux;
if(tracker) {
path.type = hit_a_receiver ? SSOL_PATH_SUCCESS : SSOL_PATH_MISSING;
}
}
+ /* Now that the sample ends successfully, record MC weights */
+ ACCUM_WEIGHT(pt.mc_samp->cos_factor, pt.cos_factor);
+ ACCUM_WEIGHT(thread_ctx->cos_factor, pt.cos_factor);
+ #undef ACCUM_WEIGHT
+
+ /* Check conservation of energy at the realisation level */
+ ASSERT((double)depth*DBL_EPSILON*pt.initial_flux >= fabs(pt.energy_loss));
+exit:
if(tracker) {
- res = path_register_and_clear(&thread_ctx->paths, &path);
- if(res != RES_OK) goto error;
+ res_T tmp_res = path_register_and_clear(&thread_ctx->paths, &path);
+ if(tmp_res != RES_OK && res == RES_OK) {
+ res = tmp_res;
+ goto error;
+ }
}
-exit:
ssol_medium_clear(&in_medium);
ssol_medium_clear(&out_medium);
if(tracker) path_release(&path);
return res;
error:
+ if (tracker) {
+ path.type = SSOL_PATH_ERROR;
+ }
goto exit;
}
+static void
+cancel_mc_receiver_1side
+ (struct mc_receiver_1side* rcv,
+ size_t irealisation)
+{
+ mc_data_cancel(&rcv->incoming_flux, irealisation);
+ mc_data_cancel(&rcv->incoming_if_no_atm_loss, irealisation);
+ mc_data_cancel(&rcv->incoming_if_no_field_loss, irealisation);
+ mc_data_cancel(&rcv->incoming_lost_in_field, irealisation);
+ mc_data_cancel(&rcv->incoming_lost_in_atmosphere, irealisation);
+ mc_data_cancel(&rcv->absorbed_flux, irealisation);
+ mc_data_cancel(&rcv->absorbed_if_no_atm_loss, irealisation);
+ mc_data_cancel(&rcv->absorbed_if_no_field_loss, irealisation);
+ mc_data_cancel(&rcv->absorbed_lost_in_field, irealisation);
+ mc_data_cancel(&rcv->absorbed_lost_in_atmosphere, irealisation);
+}
+
+static void
+cancel_mc
+ (struct thread_context* thread_ctx,
+ size_t irealisation)
+{
+ struct htable_receiver_iterator r_it, r_end;
+ struct htable_sampled_iterator s_it, s_end;
+
+ /* Cancel global MC estimations */
+ mc_data_cancel(&thread_ctx->cos_factor, irealisation);
+ mc_data_cancel(&thread_ctx->absorbed_by_atmosphere, irealisation);
+ mc_data_cancel(&thread_ctx->absorbed_by_receivers, irealisation);
+ mc_data_cancel(&thread_ctx->other_absorbed, irealisation);
+ mc_data_cancel(&thread_ctx->missing, irealisation);
+ mc_data_cancel(&thread_ctx->shadowed, irealisation);
+
+ /* Cancel receiver MC estimations */
+ htable_receiver_begin(&thread_ctx->mc_rcvs, &r_it);
+ htable_receiver_end(&thread_ctx->mc_rcvs, &r_end);
+ while (!htable_receiver_iterator_eq(&r_it, &r_end)) {
+ struct mc_receiver* mc_rcv = htable_receiver_iterator_data_get(&r_it);
+ const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&r_it);
+
+ htable_receiver_iterator_next(&r_it);
+
+ if (inst->receiver_mask & (int)SSOL_FRONT) {
+ cancel_mc_receiver_1side(&mc_rcv->front, irealisation);
+ }
+ if (inst->receiver_mask & (int)SSOL_BACK) {
+ cancel_mc_receiver_1side(&mc_rcv->back, irealisation);
+ }
+ }
+ /* Cancel sampled instance MC estimations */
+ htable_sampled_begin(&thread_ctx->mc_samps, &s_it);
+ htable_sampled_end(&thread_ctx->mc_samps, &s_end);
+ while (!htable_sampled_iterator_eq(&s_it, &s_end)) {
+ struct mc_sampled* mc_samp = htable_sampled_iterator_data_get(&s_it);
+ htable_sampled_iterator_next(&s_it);
+
+ mc_data_cancel(&mc_samp->cos_factor, irealisation);
+ mc_data_cancel(&mc_samp->shadowed, irealisation);
+
+ /* dst->by_receiver += src->by_receiver; */
+ htable_receiver_begin(&mc_samp->mc_rcvs, &r_it);
+ htable_receiver_end(&mc_samp->mc_rcvs, &r_end);
+ while (!htable_receiver_iterator_eq(&r_it, &r_end)) {
+ struct mc_receiver* mc_rcv = htable_receiver_iterator_data_get(&r_it);
+ const struct ssol_instance* inst = *htable_receiver_iterator_key_get(&r_it);
+ htable_receiver_iterator_next(&r_it);
+
+ if (inst->receiver_mask & (int)SSOL_FRONT) {
+ cancel_mc_receiver_1side(&mc_rcv->front, irealisation);
+ }
+ if (inst->receiver_mask & (int)SSOL_BACK) {
+ cancel_mc_receiver_1side(&mc_rcv->back, irealisation);
+ }
+ }
+ }
+}
+
/*******************************************************************************
* Exported functions
******************************************************************************/
@@ -1005,7 +1097,7 @@ ssol_solve
res = scene_check(scn, FUNC_NAME);
if(res != RES_OK) goto error;
-
+
/* init air properties */
if(scn->atmosphere)
ssol_data_copy(&scn->air.absorption, &scn->atmosphere->absorption);
@@ -1065,7 +1157,10 @@ ssol_solve
/* Execute a MC experiment */
res_local = trace_radiative_path((size_t)i, sampled_area_proxy, thread_ctx,
scn, view_samp, view_rt, ran_sun_dir, ran_sun_wl, path_tracker);
-
+ if(res_local != RES_OK) {
+ /* Cancel partial MC results */
+ cancel_mc(thread_ctx, (size_t)i);
+ }
if(res_local == RES_BAD_OP) {
if(ATOMIC_INCR(&nfailures) >= max_failures) {
log_error(scn->dev, "Too many unexpected radiative paths.\n");
diff --git a/src/test_ssol_solver4.c b/src/test_ssol_solver4.c
@@ -154,11 +154,9 @@ main(int argc, char** argv)
m2 = m1;
std2 = std1;
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
- printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE);
- printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE);
- printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE);
+ PRINT_GLOBAL(mc_global);
CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1);
- CHECK(eq_eps(mc_global.missing.E, 400 * DNI_cos, 1e-2), 1); /* nothing absorbed */
+ CHECK(eq_eps(mc_global.missing.E, 400 * DNI_cos, 1e-2), 1); /* virtual => 100 % missing */
CHECK(GET_MC_RCV(estimator, target1, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target1) = %g +/- %g\n",
mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
diff --git a/src/test_ssol_solver9.c b/src/test_ssol_solver9.c
@@ -141,9 +141,8 @@ main(int argc, char** argv)
#define DNI_TGT_S (DNI * TGT_X * TGT_Y)
#define DNI_S (DNI * SZ * SZ)
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
- printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE);
- printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE);
- printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE);
+ PRINT_GLOBAL(mc_global);
+ CHECK(eq_eps(mc_global.cos_factor.E, 1./3., 3 * mc_global.cos_factor.SE), 1);
CHECK(eq_eps(mc_global.shadowed.E, DNI_S, 3 * mc_global.shadowed.SE), 1);
CHECK(eq_eps(mc_global.missing.E, MMAX(DNI_S, DNI_TGT_S),
3 * mc_global.missing.SE), 1);
