commit 443e383981221ca0213d60caedd46123a0336138
parent 3cfa90f7822b98c68b202d5c72363da08623e3c5
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 3 May 2021 16:27:35 +0200
Add the compute_4v_s function
Diffstat:
3 files changed, 187 insertions(+), 2 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -54,6 +54,7 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SGS_FILES_SRC
sgs_args.c
sgs.c
+ sgs_compute_4v_s.c
sgs_geometry.c
sgs_geometry_box.c
sgs_geometry_slope.c
diff --git a/src/sgs.c b/src/sgs.c
@@ -50,7 +50,7 @@ sgs_create
goto error;
}
sgs->allocator = allocator;
- sgs->nrealisations = (size_t)sgs->nrealisations;
+ sgs->nrealisations = (size_t)args->nrealisations;
sgs->nthreads = args->nthreads;
sgs->verbose = args->verbose;
sgs->dump_geometry = args->dump_geometry;
@@ -144,7 +144,8 @@ sgs_run(struct sgs* sgs)
res = sgs_geometry_dump_vtk(sgs->geom, stdout);
if(res != RES_OK) goto error;
} else {
- /* TODO Monte Carlo computation */
+ res = compute_4v_s(sgs);
+ if(res != RES_OK) goto error;
}
exit:
diff --git a/src/sgs_compute_4v_s.c b/src/sgs_compute_4v_s.c
@@ -0,0 +1,183 @@
+/* Copyright (C) 2021
+ * CNRS/RAPSODEE,
+ * CNRS/LMAP,
+ * |Meso|Star> (contact@meso-star.com),
+ * UPS/Laplace.
+ *
+ * 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 "sgs_c.h"
+#include "sgs_geometry.h"
+#include "sgs_log.h"
+
+#include <star/smc.h>
+#include <star/ssp.h>
+
+#include <rsys/double2.h>
+#include <rsys/double3.h>
+#include <rsys/cstr.h>
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+/* Reflect the vector V wrt the normal N. By convention V points outward the
+ * surface. */
+static INLINE double*
+reflect(double res[3], const double V[3], const double N[3])
+{
+ double tmp[3];
+ double cos_V_N;
+ ASSERT(res && V && N);
+ ASSERT(d3_is_normalized(V) && d3_is_normalized(N));
+ cos_V_N = d3_dot(V, N);
+ d3_muld(tmp, N, 2*cos_V_N);
+ d3_sub(res, tmp, V);
+ return res;
+}
+
+static res_T
+realisation(void* weight, struct ssp_rng* rng, const unsigned ithread, void* ctx)
+{
+ struct sgs_fragment frag = SGS_FRAGMENT_NULL;
+ struct sgs_hit hit = SGS_HIT_NULL;
+
+ double vec[3] = {0, 0, 0};
+ double pos[3] = {0, 0, 0};
+ double dir[3] = {0, 0, 0};
+ double range[2] = {0, 0};
+
+ struct sgs_s3d_position* pos_from = NULL;
+ struct sgs* sgs = ctx;
+
+ double len = 0;
+
+ int sampling_mask = 0;
+ res_T res = RES_OK;
+
+ ASSERT(weight && rng && ctx);
+ (void)ithread;
+
+ sampling_mask = sgs_geometry_get_sampling_mask(sgs->geom);
+
+ /* Sample the emissive surface */
+ sgs_geometry_sample(sgs->geom, rng, &frag);
+ pos_from = &frag.s3d_pos;
+
+ /* Cosine weighted sampling of the hemisphere around the sampled normal */
+ ssp_ran_hemisphere_cos(rng, frag.normal, dir, NULL);
+
+ pos[0] = frag.position[0];
+ pos[1] = frag.position[1];
+ pos[2] = frag.position[2];
+
+ range[0] = 0;
+ range[1] = INF;
+
+ for(;;) {
+ sgs_geometry_trace_ray(sgs->geom, pos, dir, range, pos_from, &hit);
+ if(SGS_HIT_NONE(&hit)) { /* Unexpected error */
+ res = RES_BAD_OP;
+ goto error;
+ }
+
+ len += hit.distance;
+
+ if(((int)hit.surface_type & sampling_mask) != 0)
+ break; /* Stop random walk */
+
+ /* Move to the intersection */
+ d3_muld(vec, dir, hit.distance);
+ d3_add(pos, pos, vec);
+
+ /* Reflect the direction with respect to the intersection normal */
+ d3_minus(dir, dir); /* Ensure reflect convention */
+ reflect(dir, dir, hit.normal);
+
+ pos_from = &hit.s3d_pos;
+ }
+exit:
+ SMC_DOUBLE(weight) = len;
+ return res;
+error:
+ len = NaN;
+ goto exit;
+}
+
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+res_T
+compute_4v_s(struct sgs* sgs)
+{
+ struct smc_integrator integrator;
+ struct smc_estimator_status status;
+ struct smc_device* smc = NULL;
+ struct smc_estimator* estimator = NULL;
+ double V;
+ double S;
+ res_T res = RES_OK;
+ ASSERT(sgs);
+
+ /* Create the Star-MonteCarlo device */
+ res = smc_device_create
+ (&sgs->logger, sgs->allocator, sgs->nthreads, &ssp_rng_mt19937_64, &smc);
+ if(res != RES_OK) {
+ sgs_log_err(sgs, "Could not create the Star-MonteCarlo device -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ /* Setup the integrator */
+ integrator.integrand = realisation;
+ integrator.type = &smc_double;
+ integrator.max_realisations = sgs->nrealisations;
+ integrator.max_failures = (size_t)((double)sgs->nrealisations*0.01);
+
+ /* Run Monte-Carlo integration */
+ res = smc_solve(smc, &integrator, sgs, &estimator);
+ if(res != RES_OK) {
+ sgs_log_err(sgs, "Error while solving 4V/S -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ /* Retrieve the estimated result */
+ SMC(estimator_get_status(estimator, &status));
+ if(status.NF == integrator.max_failures) {
+ sgs_log_err(sgs, "Too many failures: %lu\n", (unsigned long)status.NF);
+ res = RES_BAD_OP;
+ goto error;
+ }
+
+ /* Fetch the reference values */
+ V = sgs_geometry_compute_volume(sgs->geom);
+ S = sgs_geometry_compute_sampling_area(sgs->geom);
+ sgs_log(sgs, "V = %g; S = %g\n", V, S);
+
+ /* Print the result */
+ sgs_log(sgs, "4V/S = %g ~ %g +/- %g; #failures = %lu/%lu\n",
+ 4*V/S,
+ SMC_DOUBLE(status.E),
+ SMC_DOUBLE(status.SE),
+ (unsigned long)status.NF,
+ (unsigned long)sgs->nrealisations);
+
+exit:
+ if(smc) SMC(device_ref_put(smc));
+ if(estimator) SMC(estimator_ref_put(estimator));
+ return res;
+error:
+ goto exit;
+}
