solstice-solver

commit 3fcaff89737ddb646f49eb440f06fad3187d4009
parent 6e99b59312815340851f3a4a7188ff6369ae6b73
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Tue, 12 Sep 2017 11:54:37 +0200

BugFix: prevent realisations in error to change MC results.

By adding a cancel mechanism on MC results.

Diffstat:
M
src/ssol_estimator_c.h
 | 
8
++++++++
M
src/ssol_solver.c
 | 
82
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
M
src/test_ssol_solver4.c
 | 
6
++----

3 files changed, 91 insertions(+), 5 deletions(-)
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_solver.c b/src/ssol_solver.c
@@ -959,6 +959,83 @@ 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
  ******************************************************************************/
@@ -1069,7 +1146,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);