commit e6e567e2a6c398c6e5b65f9fb2c1356b78949d8d
parent 922814aeb8a1977f27ee6897a5160d2c957e71da
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 9 Feb 2024 10:43:30 +0100
Merge branch 'feature_solve_probe_list' into develop
Diffstat:
13 files changed, 1121 insertions(+), 177 deletions(-)
diff --git a/Makefile b/Makefile
@@ -212,7 +212,8 @@ TEST_SRC_MPI =\
src/test_sdis_solve_boundary.c\
src/test_sdis_solve_boundary_flux.c\
src/test_sdis_solve_probe2.c\
- src/test_sdis_solve_probe_list.c
+ src/test_sdis_solve_probe_list.c\
+ src/test_sdis_solve_probe_boundary_list.c
TEST_OBJ =\
$(TEST_SRC:.c=.o)\
$(TEST_SRC_MPI:.c=.o)\
@@ -430,18 +431,21 @@ test_sdis_solve_medium_2d \
src/test_sdis_compute_power.d \
src/test_sdis_solve_camera.d \
src/test_sdis_solve_medium.d \
+src/test_sdis_solve_probe_boundary_list.d \
: config.mk sdis-local.pc
@$(CC) $(TEST_CFLAGS_MPI) $(S3DUT_CFLAGS) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
src/test_sdis_compute_power.o \
src/test_sdis_solve_camera.o \
src/test_sdis_solve_medium.o \
+src/test_sdis_solve_probe_boundary_list.o \
: config.mk sdis-local.pc
$(CC) $(TEST_CFLAGS_MPI) $(S3DUT_CFLAGS) -c $(@:.o=.c) -o $@
test_sdis_compute_power \
test_sdis_solve_camera \
test_sdis_solve_medium \
+test_sdis_solve_probe_boundary_list \
: config.mk sdis-local.pc $(LIBNAME) src/test_sdis_utils.o
$(CC) $(TEST_CFLAGS_MPI) $(S3DUT_CFLAGS) -o $@ src/$@.o $(TEST_LIBS_MPI) $(S3DUT_LIBS)
diff --git a/src/sdis.c b/src/sdis.c
@@ -476,10 +476,10 @@ compute_process_index_range
per_process_indices = nindices / (size_t)dev->mpi_nprocs;
range[0] = per_process_indices * (size_t)dev->mpi_rank;
- range[1] = range[0] + per_process_indices; /* Upper bound is _exclusive */
+ range[1] = range[0] + per_process_indices; /* Upper bound is _exclusive_ */
ASSERT(range[0] <= range[1]);
- /* Set the remaining number of indexes that are not managed by one process */
+ /* Set the remaining number of indices that are not managed by one process */
remaining_indices =
nindices - per_process_indices * (size_t)dev->mpi_nprocs;
diff --git a/src/sdis.h b/src/sdis.h
@@ -553,6 +553,27 @@ static const struct sdis_solve_probe_boundary_args
SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT =
SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT__;
+/* Input arguments of the solve function that distributes the calculations of
+ * several boundary probes rather than the realizations of a probe */
+struct sdis_solve_probe_boundary_list_args {
+ struct sdis_solve_probe_boundary_args* probes; /* List of probes to compute */
+ size_t nprobes; /* Total number of probes */
+
+ /* State/type of the RNG to use for the list of probes to calculate.
+ * The state/type defines per probe is ignored */
+ struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
+ enum ssp_rng_type rng_type; /* RNG type to use if `rng_state' is NULL */
+};
+#define SDIS_SOLVE_PROBE_BOUNDARY_LIST_ARGS_DEFAULT__ { \
+ NULL, /* List of probes */ \
+ 0, /* #probes */ \
+ NULL, /* RNG state */ \
+ SSP_RNG_THREEFRY /* RNG type */ \
+}
+static const struct sdis_solve_probe_boundary_list_args
+SDIS_SOLVE_PROBE_BOUNDARY_LIST_ARGS_DEFAULT =
+ SDIS_SOLVE_PROBE_BOUNDARY_LIST_ARGS_DEFAULT__;
+
struct sdis_solve_boundary_args {
size_t nrealisations; /* #realisations */
const size_t* primitives; /* List of boundary primitives to handle */
@@ -1424,18 +1445,6 @@ sdis_solve_probe
const struct sdis_solve_probe_args* args,
struct sdis_estimator** estimator);
-/* Calculate temperature for a list of probe points. Unlike its
- * single-probe counterpart, this function parallelizes the list of
- * probes, rather than calculating a single probe. Calling this function
- * is therefore more advantageous in terms of load distribution when the
- * number of probe points to be evaluated is large compared to the cost
- * of calculating a single probe point. */
-SDIS_API res_T
-sdis_solve_probe_list
- (struct sdis_scene* scn,
- const struct sdis_solve_probe_list_args* args,
- struct sdis_estimator_buffer** buf);
-
SDIS_API res_T
sdis_solve_probe_boundary
(struct sdis_scene* scn,
@@ -1486,6 +1495,27 @@ sdis_compute_power
struct sdis_estimator** estimator);
/*******************************************************************************
+ * Solvers of a list of probes
+ *
+ * Unlike their single-probe counterpart, this function parallelizes the list of
+ * probes, rather than calculating a single probe. Calling these functions is
+ * therefore more advantageous in terms of load distribution when the number of
+ * probes to be evaluated is large compared to the cost of calculating a single
+ * probe.
+ ******************************************************************************/
+SDIS_API res_T
+sdis_solve_probe_list
+ (struct sdis_scene* scn,
+ const struct sdis_solve_probe_list_args* args,
+ struct sdis_estimator_buffer** buf);
+
+SDIS_API res_T
+sdis_solve_probe_boundary_list
+ (struct sdis_scene* scn,
+ const struct sdis_solve_probe_boundary_list_args* args,
+ struct sdis_estimator_buffer** buf);
+
+/*******************************************************************************
* Green solvers.
*
* Note that only the interfaces/media with flux/volumic power defined during
diff --git a/src/sdis_estimator_buffer.c b/src/sdis_estimator_buffer.c
@@ -273,3 +273,8 @@ estimator_buffer_save_rng_state
return create_rng_from_rng_proxy(buf->dev, proxy, &buf->rng);
}
+/* Define the functions generic to the observable type */
+#define SDIS_X_OBS probe
+#include "sdis_estimator_buffer_X_obs.h"
+#define SDIS_X_OBS probe_boundary
+#include "sdis_estimator_buffer_X_obs.h"
diff --git a/src/sdis_estimator_buffer_X_obs.h b/src/sdis_estimator_buffer_X_obs.h
@@ -0,0 +1,118 @@
+/* Copyright (C) 2016-2023 |Méso|Star> (contact@meso-star.com)
+ *
+ * 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/>. */
+
+/* Define functions on estimator buffer generic to the observable type. */
+
+#ifndef SDIS_ESTIMATOR_BUFFER_X_OBS_H
+#define SDIS_ESTIMATOR_BUFFER_X_OBS_H
+
+#include "sdis.h"
+#include "sdis_estimator_c.h"
+#include "sdis_estimator_buffer_c.h"
+#include "sdis_misc.h"
+
+#endif /* SDIS_ESTIMATOR_BUFFER_X_OBS_H */
+
+/* Check the generic observable parameter */
+#ifndef SDIS_X_OBS
+ #error "The SDIS_X_OBS macro must be defined."
+#endif
+
+#define X_OBS(Name) CONCAT(CONCAT(Name, _), SDIS_X_OBS)
+#define SDIS_SOLVE_X_OBS_ARGS CONCAT(CONCAT(sdis_solve_, SDIS_X_OBS),_args)
+
+res_T
+X_OBS(estimator_buffer_create_from_observable_list)
+ (struct sdis_device* dev,
+ struct ssp_rng_proxy* rng_proxy,
+ const struct SDIS_SOLVE_X_OBS_ARGS obs_list_args[],
+ const struct accum* per_obs_acc_temp,
+ const struct accum* per_obs_acc_time,
+ const size_t nobs, /* #observables */
+ struct sdis_estimator_buffer** out_estim_buffer)
+{
+ /* Accumulators throughout the buffer */
+ struct accum acc_temp = ACCUM_NULL;
+ struct accum acc_time = ACCUM_NULL;
+ size_t nrealisations = 0;
+
+ struct sdis_estimator_buffer* estim_buf = NULL;
+ size_t iobs = 0;
+ res_T res = RES_OK;
+
+ ASSERT(dev && rng_proxy);
+ ASSERT(obs_list_args || !nobs);
+ ASSERT(per_obs_acc_time && per_obs_acc_time && out_estim_buffer);
+
+ res = estimator_buffer_create(dev, nobs, 1, &estim_buf);
+ if(res != RES_OK) {
+ log_err(dev, "Unable to allocate the estimator buffer.\n");
+ goto error;
+ }
+
+ FOR_EACH(iobs, 0, nobs) {
+ const struct SDIS_SOLVE_X_OBS_ARGS* obs_args = NULL;
+ const struct accum* obs_acc_temp = NULL;
+ const struct accum* obs_acc_time = NULL;
+ struct sdis_estimator* estim = NULL;
+
+ /* Get observable data */
+ obs_args = obs_list_args + iobs;
+ obs_acc_temp = per_obs_acc_temp + iobs;
+ obs_acc_time = per_obs_acc_time + iobs;
+ ASSERT(obs_acc_temp->count == obs_acc_time->count);
+
+ /* Setup the estimator of the current observable */
+ estim = estimator_buffer_grab(estim_buf, iobs, 0);
+ estimator_setup_realisations_count
+ (estim, obs_args->nrealisations, obs_acc_temp->count);
+ estimator_setup_temperature
+ (estim, obs_acc_temp->sum, obs_acc_temp->sum2);
+ estimator_setup_realisation_time
+ (estim, obs_acc_time->sum, obs_acc_time->sum2);
+
+ /* Update global accumulators */
+ acc_temp.sum += obs_acc_temp->sum;
+ acc_temp.sum2 += obs_acc_temp->sum2;
+ acc_temp.count += obs_acc_temp->count;
+ acc_time.sum += obs_acc_time->sum;
+ acc_time.sum2 += obs_acc_time->sum2;
+ acc_time.count += obs_acc_time->count;
+ nrealisations += obs_args->nrealisations;
+ }
+
+ ASSERT(acc_temp.count == acc_time.count);
+
+ /* Setup global estimator */
+ estimator_buffer_setup_realisations_count
+ (estim_buf, nrealisations, acc_temp.count);
+ estimator_buffer_setup_temperature
+ (estim_buf, acc_temp.sum, acc_temp.sum2);
+ estimator_buffer_setup_realisation_time
+ (estim_buf, acc_time.sum, acc_time.sum2);
+
+ res = estimator_buffer_save_rng_state(estim_buf, rng_proxy);
+ if(res != RES_OK) goto error;
+
+exit:
+ *out_estim_buffer = estim_buf;
+ return res;
+error:
+ goto exit;
+}
+
+#undef X_OBS
+#undef SDIS_SOLVE_X_OBS_ARGS
+#undef SDIS_X_OBS
diff --git a/src/sdis_estimator_buffer_c.h b/src/sdis_estimator_buffer_c.h
@@ -59,5 +59,25 @@ estimator_buffer_save_rng_state
(struct sdis_estimator_buffer* buf,
const struct ssp_rng_proxy* proxy);
+extern LOCAL_SYM res_T
+estimator_buffer_create_from_observable_list_probe
+ (struct sdis_device* dev,
+ struct ssp_rng_proxy* rng_proxy,
+ const struct sdis_solve_probe_args obs_list_args[],
+ const struct accum* per_obs_acc_temp,
+ const struct accum* per_obs_acc_time,
+ const size_t nobs, /* #observables */
+ struct sdis_estimator_buffer** out_estim_buffer);
+
+extern LOCAL_SYM res_T
+estimator_buffer_create_from_observable_list_probe_boundary
+ (struct sdis_device* dev,
+ struct ssp_rng_proxy* rng_proxy,
+ const struct sdis_solve_probe_boundary_args obs_list_args[],
+ const struct accum* per_obs_acc_temp,
+ const struct accum* per_obs_acc_time,
+ const size_t nobs, /* #observables */
+ struct sdis_estimator_buffer** out_estim_buffer);
+
#endif /* SDIS_ESTIMATOR_BUFFER_C_H */
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -59,7 +59,7 @@ sdis_solve_probe
res_T
sdis_solve_probe_list
(struct sdis_scene* scn,
- const struct sdis_solve_probe_list_args args[],
+ const struct sdis_solve_probe_list_args* args,
struct sdis_estimator_buffer** out_buf)
{
if(!scn) return RES_BAD_ARG;
@@ -99,6 +99,20 @@ sdis_solve_probe_boundary
}
res_T
+sdis_solve_probe_boundary_list
+ (struct sdis_scene* scn,
+ const struct sdis_solve_probe_boundary_list_args* args,
+ struct sdis_estimator_buffer** out_buf)
+{
+ if(!scn) return RES_BAD_ARG;
+ if(scene_is_2d(scn)) {
+ return solve_probe_boundary_list_2d(scn, args, out_buf);
+ } else {
+ return solve_probe_boundary_list_3d(scn, args, out_buf);
+ }
+}
+
+res_T
sdis_solve_probe_boundary_green_function
(struct sdis_scene* scn,
const struct sdis_solve_probe_boundary_args* args,
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -101,84 +101,6 @@ check_solve_probe_list_args
return RES_OK;
}
-static res_T
-setup_estimator_buffer
- (struct sdis_device* dev,
- struct ssp_rng_proxy* rng_proxy,
- const struct sdis_solve_probe_list_args* solve_args,
- const struct accum* per_probe_acc_temp,
- const struct accum* per_probe_acc_time,
- struct sdis_estimator_buffer** out_estim_buffer)
-{
- /* Accumulators throughout the buffer */
- struct accum acc_temp = ACCUM_NULL;
- struct accum acc_time = ACCUM_NULL;
- size_t nrealisations = 0;
-
- struct sdis_estimator_buffer* estim_buf = NULL;
- size_t iprobe = 0;
- res_T res = RES_OK;
-
- ASSERT(dev && rng_proxy && solve_args);
- ASSERT(per_probe_acc_time && per_probe_acc_time && out_estim_buffer);
-
- res = estimator_buffer_create(dev, solve_args->nprobes, 1, &estim_buf);
- if(res != RES_OK) {
- log_err(dev, "Unable to allocate the estimator buffer.\n");
- goto error;
- }
-
- FOR_EACH(iprobe, 0, solve_args->nprobes) {
- const struct sdis_solve_probe_args* probe = NULL;
- const struct accum* probe_acc_temp = NULL;
- const struct accum* probe_acc_time = NULL;
- struct sdis_estimator* estim = NULL;
-
- /* Get probe data */
- probe = solve_args->probes + iprobe;
- probe_acc_temp = per_probe_acc_temp + iprobe;
- probe_acc_time = per_probe_acc_time + iprobe;
- ASSERT(probe_acc_temp->count == probe_acc_time->count);
-
- /* Setup probe estimator */
- estim = estimator_buffer_grab(estim_buf, iprobe, 0);
- estimator_setup_realisations_count
- (estim, probe->nrealisations, probe_acc_temp->count);
- estimator_setup_temperature
- (estim, probe_acc_temp->sum, probe_acc_temp->sum2);
- estimator_setup_realisation_time
- (estim, probe_acc_time->sum, probe_acc_time->sum2);
-
- /* Update global accumulators */
- acc_temp.sum += probe_acc_temp->sum;
- acc_temp.sum2 += probe_acc_temp->sum2;
- acc_temp.count += probe_acc_temp->count;
- acc_time.sum += probe_acc_time->sum;
- acc_time.sum2 += probe_acc_time->sum2;
- acc_time.count += probe_acc_time->count;
- nrealisations += probe->nrealisations;
- }
-
- ASSERT(acc_temp.count == acc_time.count);
-
- /* Setup global estimator */
- estimator_buffer_setup_realisations_count
- (estim_buf, nrealisations, acc_temp.count);
- estimator_buffer_setup_temperature
- (estim_buf, acc_temp.sum, acc_temp.sum2);
- estimator_buffer_setup_realisation_time
- (estim_buf, acc_time.sum, acc_time.sum2);
-
- res = estimator_buffer_save_rng_state(estim_buf, rng_proxy);
- if(res != RES_OK) goto error;
-
-exit:
- *out_estim_buffer = estim_buf;
- return res;
-error:
- goto exit;
-}
-
#endif /* SDIS_SOLVE_PROBE_XD_H */
static res_T
@@ -626,9 +548,6 @@ XD(solve_probe_list)
goto post_sync;
}
- /* Begin time registration of the computation */
- time_current(&time0);
-
/* Allocate the list of accumulators per probe. On the master process,
* allocate a complete list in which the accumulators of all processes will be
* stored. */
@@ -641,6 +560,9 @@ XD(solve_probe_list)
goto error;
}
+ /* Begin time registration of the computation */
+ time_current(&time0);
+
/* Here we go! Calculation of probe list */
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
@@ -726,8 +648,9 @@ post_sync:
log_info(scn->dev, "Probes accumulator gathered in %s.\n", buf);
if(is_master_process) {
- res = setup_estimator_buffer(scn->dev, rng_proxy, args, per_probe_acc_temp,
- per_probe_acc_time, &estim_buf);
+ res = estimator_buffer_create_from_observable_list_probe
+ (scn->dev, rng_proxy, args->probes, per_probe_acc_temp,
+ per_probe_acc_time, args->nprobes, &estim_buf);
if(res != RES_OK) goto error;
}
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -75,6 +75,40 @@ check_solve_probe_boundary_args
}
static INLINE res_T
+check_solve_probe_boundary_list_args
+ (struct sdis_device* dev,
+ const struct sdis_solve_probe_boundary_list_args* args)
+{
+ size_t iprobe = 0;
+
+ if(!args) return RES_BAD_ARG;
+
+ /* Check the list of probes */
+ if(!args->probes || !args->nprobes) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check the RNG type */
+ if(!args->rng_state && args->rng_type >= SSP_RNG_TYPES_COUNT__) {
+ return RES_BAD_ARG;
+ }
+
+ FOR_EACH(iprobe, 0, args->nprobes) {
+ const res_T res = check_solve_probe_boundary_args(args->probes+iprobe);
+ if(res != RES_OK) return res;
+
+ if(args->probes[iprobe].register_paths != SDIS_HEAT_PATH_NONE) {
+ log_warn(dev,
+ "Unable to save paths for probe boundary %lu. "
+ "Saving path is not supported when solving multiple probes\n",
+ (unsigned long)iprobe);
+ }
+ }
+
+ return RES_OK;
+}
+
+static INLINE res_T
check_solve_probe_boundary_flux_args
(const struct sdis_solve_probe_boundary_flux_args* args)
{
@@ -109,6 +143,66 @@ check_solve_probe_boundary_flux_args
#endif /* SDIS_SOLVE_PROBE_BOUNDARY_XD_H */
+static res_T
+XD(solve_one_probe_boundary)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_solve_probe_boundary_args* args,
+ struct accum* acc_temp,
+ struct accum* acc_time)
+{
+ size_t irealisation = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && rng && check_solve_probe_boundary_args(args) == RES_OK);
+
+ *acc_temp = ACCUM_NULL;
+ *acc_time = ACCUM_NULL;
+
+ FOR_EACH(irealisation, 0, args->nrealisations) {
+ struct boundary_realisation_args realis_args = BOUNDARY_REALISATION_ARGS_NULL;
+ double w = NaN; /* MC weight */
+ double usec = 0; /* Time of a realisation */
+ double time = 0; /* Sampled observation time */
+ struct time t0, t1; /* Register the time spent solving a realisation */
+
+ /* Begin time registration of the realisation */
+ time_current(&t0);
+
+ /* Sample observation time */
+ time = sample_time(rng, args->time_range);
+
+ /* Run a realisation */
+ realis_args.rng = rng;
+ realis_args.iprim = args->iprim;
+ realis_args.time = time;
+ realis_args.picard_order = args->picard_order;
+ realis_args.side = args->side;
+ realis_args.irealisation = irealisation;
+ realis_args.uv[0] = args->uv[0];
+#if SDIS_XD_DIMENSION == 3
+ realis_args.uv[1] = args->uv[1];
+#endif
+ res = XD(boundary_realisation)(scn, &realis_args, &w);
+ if(res != RES_OK) goto error;
+
+ /* Stop time registration */
+ time_sub(&t0, time_current(&t1), &t0);
+ usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+
+ /* Update MC weights */
+ acc_temp->sum += w;
+ acc_temp->sum2 += w*w;
+ acc_temp->count += 1;
+ acc_time->sum += usec;
+ acc_time->sum2 += usec*usec;
+ acc_time->count += 1;
+ }
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -411,6 +505,207 @@ error:
}
static res_T
+XD(solve_probe_boundary_list)
+ (struct sdis_scene* scn,
+ const struct sdis_solve_probe_boundary_list_args* args,
+ struct sdis_estimator_buffer** out_estim_buf)
+{
+ /* Time registration */
+ struct time time0, time1;
+ char buf[128]; /* Temporary buffer used to store formated time */
+
+ /* Device variable */
+ struct mem_allocator* allocator = NULL;
+
+ /* Stardis variables */
+ struct sdis_estimator_buffer* estim_buf = NULL;
+
+ /* Random Number generator */
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** per_thread_rng = NULL;
+
+ /* Probe variables */
+ size_t process_probes[2] = {0, 0}; /* Probes range managed by the process */
+ size_t process_nprobes = 0; /* Number of probes managed by the process */
+ size_t nprobes = 0;
+ struct accum* per_probe_acc_temp = NULL;
+ struct accum* per_probe_acc_time = NULL;
+
+ /* Miscellaneous */
+ int32_t* progress = NULL; /* Per process progress bar */
+ int is_master_process = 0;
+ int pcent_progress = 1; /* Percentage requiring progress update */
+ int64_t i = 0;
+ ATOMIC nsolved_probes = 0;
+ ATOMIC res = RES_OK;
+
+ if(!scn || !out_estim_buf) { res = RES_BAD_ARG; goto error; }
+ res = check_solve_probe_boundary_list_args(scn->dev, args);
+ if(res != RES_OK) goto error;
+ res = XD(scene_check_dimensionality)(scn);
+ if(res != RES_OK) goto error;
+
+#ifdef SDIS_ENABLE_MPI
+ is_master_process = !scn->dev->use_mpi || scn->dev->mpi_rank == 0;
+#endif
+
+ allocator = scn->dev->allocator;
+
+ /* Update the progress bar every percent if escape sequences are allowed in
+ * log messages or only every 10 percent when only plain text is allowed.
+ * This reduces the number of lines of plain text printed */
+ pcent_progress = scn->dev->no_escape_sequence ? 10 : 1;
+
+ /* Create the per threads RNGs */
+ res = create_per_thread_rng
+ (scn->dev, args->rng_state, args->rng_type, &rng_proxy, &per_thread_rng);
+ if(res != RES_OK) goto error;
+
+ /* Allocate the per process progress status */
+ res = alloc_process_progress(scn->dev, &progress);
+ if(res != RES_OK) goto error;
+
+ /* Synchronise the processes */
+ process_barrier(scn->dev);
+
+ /* Define the range of probes to manage in this process */
+ process_nprobes = compute_process_index_range
+ (scn->dev, args->nprobes, process_probes);
+
+ #define PROGRESS_MSG "Solving surface probes: "
+ print_progress(scn->dev, progress, PROGRESS_MSG);
+
+ /* If there is no work to be done on this process (i.e. no probe to
+ * calculate), simply print its completion and go straight to the
+ * synchronization barrier.*/
+ if(process_nprobes == 0) {
+ progress[0] = 100;
+ print_progress_update(scn->dev, progress, PROGRESS_MSG);
+ goto post_sync;
+ }
+
+ /* Allocate the list of accumulators per probe. On the master process,
+ * allocate a complete list in which the accumulators of all processes will be
+ * stored. */
+ nprobes = is_master_process ? args->nprobes : process_nprobes;
+ per_probe_acc_temp = MEM_CALLOC(allocator, nprobes, sizeof(*per_probe_acc_temp));
+ per_probe_acc_time = MEM_CALLOC(allocator, nprobes, sizeof(*per_probe_acc_time));
+ if(!per_probe_acc_temp || !per_probe_acc_time) {
+ log_err(scn->dev, "Unable to allocate the list of accumulators per probe.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ /* Begin time registration of the computation */
+ time_current(&time0);
+
+ /* Calculation of probe list */
+ #pragma omp parallel for schedule(static)
+ for(i = 0; i < (int64_t)process_nprobes; ++i) {
+ /* Thread */
+ const int ithread = omp_get_thread_num(); /* Thread ID */
+ struct ssp_rng* rng = per_thread_rng[ithread]; /* RNG of the thread */
+
+ /* Probe */
+ struct accum* probe_acc_temp = NULL;
+ struct accum* probe_acc_time = NULL;
+ const struct sdis_solve_probe_boundary_args* probe_args = NULL;
+ const size_t iprobe = process_probes[0] + (size_t)i; /* Probe ID */
+
+ /* Miscellaneous */
+ size_t n = 0; /* Number of solved probes */
+ int pcent = 0; /* Current progress */
+ res_T res_local = RES_OK;
+
+ if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+
+ /* Retrieve the probe arguments */
+ probe_args = &args->probes[iprobe];
+
+ /* Retrieve the probe accumulators */
+ probe_acc_temp = &per_probe_acc_temp[i];
+ probe_acc_time = &per_probe_acc_time[i];
+
+ res_local = XD(solve_one_probe_boundary)
+ (scn, rng, probe_args, probe_acc_temp, probe_acc_time);
+ if(res_local != RES_OK) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+
+ /* Update progress */
+ n = (size_t)ATOMIC_INCR(&nsolved_probes);
+ pcent = (int)((double)n * 100.0 / (double)process_nprobes + 0.5/*round*/);
+
+ #pragma omp critical
+ if(pcent/pcent_progress > progress[0]/pcent_progress) {
+ progress[0] = pcent;
+ print_progress_update(scn->dev, progress, PROGRESS_MSG);
+ }
+ }
+
+post_sync:
+ /* Synchronise processes */
+ process_barrier(scn->dev);
+
+ res = gather_res_T(scn->dev, (res_T)res);
+ if(res != RES_OK) goto error;
+
+ print_progress_completion(scn->dev, progress, PROGRESS_MSG);
+ #undef PROGRESS_MSG
+
+ /* Report computatio time */
+ time_sub(&time0, time_current(&time1), &time0);
+ time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf));
+ log_info(scn->dev, "%lu surface probes solved in %s.\n",
+ (unsigned long)args->nprobes, buf);
+
+ /* Gather the RNG proxy sequence IDs and ensure that the RNG proxy state of
+ * the master process is greater than the RNG proxy state of all other
+ * processes */
+ res = gather_rng_proxy_sequence_id(scn->dev, rng_proxy);
+ if(res != RES_OK) goto error;
+
+ time_current(&time0);
+
+ /* Gather the list of accumulators */
+ #define GATHER_ACCUMS_LIST(Msg, Acc) { \
+ res = gather_accumulators_list \
+ (scn->dev, Msg, args->nprobes, process_probes, per_probe_acc_##Acc); \
+ if(res != RES_OK) goto error; \
+ } (void)0
+ GATHER_ACCUMS_LIST(MPI_SDIS_MSG_ACCUM_TEMP, temp);
+ GATHER_ACCUMS_LIST(MPI_SDIS_MSG_ACCUM_TIME, time);
+ #undef GATHER_ACCUMS_LIST
+
+ time_sub(&time0, time_current(&time1), &time0);
+ time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf));
+ log_info(scn->dev, "Probes accumulator gathered in %s.\n", buf);
+
+ if(is_master_process) {
+ res = estimator_buffer_create_from_observable_list_probe_boundary
+ (scn->dev, rng_proxy, args->probes, per_probe_acc_temp,
+ per_probe_acc_time, args->nprobes, &estim_buf);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ if(per_thread_rng) release_per_thread_rng(scn->dev, per_thread_rng);
+ if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ if(per_probe_acc_temp) MEM_RM(allocator, per_probe_acc_temp);
+ if(per_probe_acc_time) MEM_RM(allocator, per_probe_acc_time);
+ if(progress) free_process_progress(scn->dev, progress);
+ if(out_estim_buf) *out_estim_buf = estim_buf;
+ return (res_T)res;
+error:
+ if(estim_buf) {
+ SDIS(estimator_buffer_ref_put(estim_buf));
+ estim_buf = NULL;
+ }
+ goto exit;
+}
+
+static res_T
XD(solve_probe_boundary_flux)
(struct sdis_scene* scn,
const struct sdis_solve_probe_boundary_flux_args* args,
diff --git a/src/test_sdis_draw_external_flux.c b/src/test_sdis_draw_external_flux.c
@@ -14,10 +14,10 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include "sdis.h"
+#include "test_sdis_mesh.h"
#include "test_sdis_utils.h"
#include <star/s3dut.h>
-#include <rsys/stretchy_array.h>
#define IMG_WIDTH 320
#define IMG_HEIGHT 240
@@ -37,80 +37,6 @@
/*******************************************************************************
* Geometries
******************************************************************************/
-struct mesh {
- double* positions; /* List of double 3 */
- size_t* indices; /* List of size_t 3 */
-};
-#define MESH_NULL__ {NULL, NULL}
-static const struct mesh MESH_NULL = MESH_NULL__;
-
-static INLINE void
-mesh_init(struct mesh* mesh)
-{
- CHK(mesh);
- *mesh = MESH_NULL;
-}
-
-static INLINE void
-mesh_release(struct mesh* mesh)
-{
- CHK(mesh);
- sa_release(mesh->positions);
- sa_release(mesh->indices);
-}
-
-/* Number of vertices */
-static INLINE size_t
-mesh_nvertices(const struct mesh* mesh)
-{
- CHK(mesh);
- return sa_size(mesh->positions) / 3/* #coords per vertex */;
-}
-
-/* Number of triangles */
-static INLINE size_t
-mesh_ntriangles(const struct mesh* mesh)
-{
- CHK(mesh);
- return sa_size(mesh->indices) / 3/* #indices per triangle */;
-}
-
-static void
-mesh_append
- (struct mesh* mesh,
- const struct s3dut_mesh_data* mesh_data,
- const double in_translate[3]) /* May be NULL */
-{
- double translate[3] = {0, 0, 0};
- double* positions = NULL;
- size_t* indices = NULL;
- size_t ivert = 0;
- size_t i = 0;
- CHK(mesh != NULL);
-
- ivert = mesh_nvertices(mesh);
- positions = sa_add(mesh->positions, mesh_data->nvertices*3);
- indices = sa_add(mesh->indices, mesh_data->nprimitives*3);
-
- if(in_translate) {
- translate[0] = in_translate[0];
- translate[1] = in_translate[1];
- translate[2] = in_translate[2];
- }
-
- FOR_EACH(i, 0, mesh_data->nvertices) {
- positions[i*3 + 0] = mesh_data->positions[i*3 + 0] + translate[0];
- positions[i*3 + 1] = mesh_data->positions[i*3 + 1] + translate[1];
- positions[i*3 + 2] = mesh_data->positions[i*3 + 2] + translate[2];
- }
-
- FOR_EACH(i, 0, mesh_data->nprimitives) {
- indices[i*3 + 0] = mesh_data->indices[i*3 + 0] + ivert;
- indices[i*3 + 1] = mesh_data->indices[i*3 + 1] + ivert;
- indices[i*3 + 2] = mesh_data->indices[i*3 + 2] + ivert;
- }
-}
-
static void
mesh_add_super_shape(struct mesh* mesh)
{
@@ -125,7 +51,8 @@ mesh_add_super_shape(struct mesh* mesh)
f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7;
OK(s3dut_create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &sshape));
OK(s3dut_mesh_get_data(sshape, &sshape_data));
- mesh_append(mesh, &sshape_data, NULL);
+ mesh_append(mesh, sshape_data.positions, sshape_data.nvertices,
+ sshape_data.indices, sshape_data.nprimitives, NULL);
OK(s3dut_mesh_ref_put(sshape));
}
@@ -140,7 +67,8 @@ mesh_add_ground(struct mesh* mesh)
OK(s3dut_create_cuboid(NULL, 20, 20, DEPTH, &ground));
OK(s3dut_mesh_get_data(ground, &ground_data));
- mesh_append(mesh, &ground_data, translate);
+ mesh_append(mesh, ground_data.positions, ground_data.nvertices,
+ ground_data.indices, ground_data.nprimitives, translate);
OK(s3dut_mesh_ref_put(ground));
#undef DEPTH
}
diff --git a/src/test_sdis_mesh.h b/src/test_sdis_mesh.h
@@ -0,0 +1,120 @@
+/* Copyright (C) 2016-2023 |Méso|Star> (contact@meso-star.com)
+ *
+ * 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/>. */
+
+#ifndef TEST_SDIS_MESH_H
+#define TEST_SDIS_MESH_H
+
+#include <rsys/rsys.h>
+#include <rsys/stretchy_array.h>
+
+struct mesh {
+ double* positions; /* List of double 3 */
+ size_t* indices; /* List of size_t 3 */
+};
+#define MESH_NULL__ {NULL, NULL}
+static const struct mesh MESH_NULL = MESH_NULL__;
+
+static INLINE void
+mesh_init(struct mesh* mesh)
+{
+ CHK(mesh);
+ *mesh = MESH_NULL;
+}
+
+static INLINE void
+mesh_release(struct mesh* mesh)
+{
+ CHK(mesh);
+ sa_release(mesh->positions);
+ sa_release(mesh->indices);
+}
+
+/* Number of vertices */
+static INLINE size_t
+mesh_nvertices(const struct mesh* mesh)
+{
+ CHK(mesh);
+ return sa_size(mesh->positions) / 3/* #coords per vertex */;
+}
+
+/* Number of triangles */
+static INLINE size_t
+mesh_ntriangles(const struct mesh* mesh)
+{
+ CHK(mesh);
+ return sa_size(mesh->indices) / 3/* #indices per triangle */;
+}
+
+static INLINE void
+mesh_append
+ (struct mesh* mesh,
+ const double* in_positions,
+ const size_t in_nvertices,
+ const size_t* in_indices,
+ const size_t in_ntriangles,
+ const double in_translate[3]) /* May be NULL */
+{
+ double translate[3] = {0, 0, 0};
+ double* positions = NULL;
+ size_t* indices = NULL;
+ size_t ivert = 0;
+ size_t i = 0;
+ CHK(mesh != NULL);
+
+ ivert = mesh_nvertices(mesh);
+ positions = sa_add(mesh->positions, in_nvertices*3);
+ indices = sa_add(mesh->indices, in_ntriangles*3);
+
+ if(in_translate) {
+ translate[0] = in_translate[0];
+ translate[1] = in_translate[1];
+ translate[2] = in_translate[2];
+ }
+
+ FOR_EACH(i, 0, in_nvertices) {
+ positions[i*3 + 0] = in_positions[i*3 + 0] + translate[0];
+ positions[i*3 + 1] = in_positions[i*3 + 1] + translate[1];
+ positions[i*3 + 2] = in_positions[i*3 + 2] + translate[2];
+ }
+
+ FOR_EACH(i, 0, in_ntriangles) {
+ indices[i*3 + 0] = in_indices[i*3 + 0] + ivert;
+ indices[i*3 + 1] = in_indices[i*3 + 1] + ivert;
+ indices[i*3 + 2] = in_indices[i*3 + 2] + ivert;
+ }
+}
+
+static INLINE void
+mesh_dump(const struct mesh* mesh, FILE* stream)
+{
+ size_t i, n;
+ CHK(mesh != NULL);
+
+ n = mesh_nvertices(mesh);
+ FOR_EACH(i, 0, n) {
+ fprintf(stream, "v %g %g %g\n", SPLIT3(mesh->positions+i*3));
+ }
+
+ n = mesh_ntriangles(mesh);
+ FOR_EACH(i, 0, n) {
+ fprintf(stream, "f %lu %lu %lu\n",
+ (unsigned long)(mesh->indices[i*3+0] + 1),
+ (unsigned long)(mesh->indices[i*3+1] + 1),
+ (unsigned long)(mesh->indices[i*3+2] + 1));
+ }
+ fflush(stream);
+}
+
+#endif /* TEST_SDIS_MESH_H */
diff --git a/src/test_sdis_solve_probe_boundary_list.c b/src/test_sdis_solve_probe_boundary_list.c
@@ -0,0 +1,488 @@
+/* Copyright (C) 2016-2023 |Méso|Star> (contact@meso-star.com)
+ *
+ * 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 "sdis.h"
+
+#include "test_sdis_utils.h"
+#include "test_sdis_mesh.h"
+
+#include <star/s3dut.h>
+#include <rsys/math.h>
+
+#include <math.h>
+
+/*
+ * The system is a trilinear profile of the temperature at steady state, i.e. at
+ * each point of the system we can calculate the temperature analytically. Two
+ * forms are immersed in this temperature field: a super shape and a sphere
+ * included in the super shape. On the Monte Carlo side, the temperature is
+ * unknown everywhere except on the surface of the super shape whose
+ * temperature is defined from the aformentionned trilinear profile. We will
+ * estimate the temperature at the sphere boundary at several probe points. We
+ * should find by Monte Carlo the temperature of the trilinear profile at the
+ * position of the probe. It's the test.
+ * /\ <-- T(x,y,z)
+ * ___/ \___
+ * T(z) \ __ /
+ * | T(y) \ / \ /
+ * |/ T=? *__/ \
+ * o--- T(x) /_ __ _\
+ * \/ \/
+ */
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static double
+trilinear_profile(const double pos[3])
+{
+ /* Range in X, Y and Z in which the trilinear profile is defined */
+ const double lower = -4;
+ const double upper = +4;
+
+ /* Upper temperature limit in X, Y and Z [K]. Lower temperature limit is
+ * implicitly 0 */
+ const double a = 333; /* Upper temperature limit in X [K] */
+ const double b = 432; /* Upper temperature limit in Y [K] */
+ const double c = 579; /* Upper temperature limit in Z [K] */
+
+ double x, y, z;
+
+ /* Check pre-conditions */
+ CHK(pos);
+ CHK(lower <= pos[0] && pos[0] <= upper);
+ CHK(lower <= pos[1] && pos[1] <= upper);
+ CHK(lower <= pos[2] && pos[2] <= upper);
+
+ x = (pos[0] - lower) / (upper - lower);
+ y = (pos[1] - lower) / (upper - lower);
+ z = (pos[2] - lower) / (upper - lower);
+ return a*x + b*y + c*z;
+}
+
+static INLINE float
+rand_canonic(void)
+{
+ return (float)rand() / (float)((unsigned)RAND_MAX+1);
+}
+
+static void
+sample_sphere(double pos[3])
+{
+ const double phi = rand_canonic() * 2 * PI;
+ const double v = rand_canonic();
+ const double cos_theta = 1 - 2 * v;
+ const double sin_theta = 2 * sqrt(v * (1 - v));
+ pos[0] = cos(phi) * sin_theta;
+ pos[1] = sin(phi) * sin_theta;
+ pos[2] = cos_theta;
+}
+
+static INLINE void
+check_intersection
+ (const double val0,
+ const double eps0,
+ const double val1,
+ const double eps1)
+{
+ double interval0[2], interval1[2];
+ double intersection[2];
+ interval0[0] = val0 - eps0;
+ interval0[1] = val0 + eps0;
+ interval1[0] = val1 - eps1;
+ interval1[1] = val1 + eps1;
+ intersection[0] = MMAX(interval0[0], interval1[0]);
+ intersection[1] = MMIN(interval0[1], interval1[1]);
+ CHK(intersection[0] <= intersection[1]);
+}
+
+static void
+mesh_add_super_shape(struct mesh* mesh)
+{
+ struct s3dut_mesh* sshape = NULL;
+ struct s3dut_mesh_data sshape_data;
+ struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL;
+ struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL;
+ const double radius = 2;
+ const unsigned nslices = 256;
+
+ f0.A = 1.5; f0.B = 1; f0.M = 11.0; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2.0;
+ f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7;
+ OK(s3dut_create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &sshape));
+ OK(s3dut_mesh_get_data(sshape, &sshape_data));
+ mesh_append(mesh, sshape_data.positions, sshape_data.nvertices,
+ sshape_data.indices, sshape_data.nprimitives, NULL);
+ OK(s3dut_mesh_ref_put(sshape));
+}
+
+static void
+mesh_add_sphere(struct mesh* mesh)
+{
+ struct s3dut_mesh* sphere = NULL;
+ struct s3dut_mesh_data sphere_data;
+ const double radius = 1;
+ const unsigned nslices = 128;
+
+ OK(s3dut_create_sphere(NULL, radius, nslices, nslices/2, &sphere));
+ OK(s3dut_mesh_get_data(sphere, &sphere_data));
+ mesh_append(mesh, sphere_data.positions, sphere_data.nvertices,
+ sphere_data.indices, sphere_data.nprimitives, NULL);
+ OK(s3dut_mesh_ref_put(sphere));
+}
+
+/*******************************************************************************
+ * Solid, i.e. medium of super shape and sphere
+ ******************************************************************************/
+#define SOLID_PROP(Prop, Val) \
+ static double \
+ solid_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+SOLID_PROP(calorific_capacity, 500.0) /* [J/K/Kg] */
+SOLID_PROP(thermal_conductivity, 25.0) /* [W/m/K] */
+SOLID_PROP(volumic_mass, 7500.0) /* [kg/m^3] */
+SOLID_PROP(temperature, -1/*<=> unknown*/) /* [K] */
+SOLID_PROP(delta, 1.0/20.0) /* [m] */
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+/*******************************************************************************
+ * Dummy environment, i.e. environment surrounding the super shape. It is
+ * defined only for Stardis compliance: in Stardis, an interface must divide 2
+ * media.
+ ******************************************************************************/
+static struct sdis_medium*
+create_dummy(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* dummy = NULL;
+
+ shader.calorific_capacity = dummy_medium_getter;
+ shader.volumic_mass = dummy_medium_getter;
+ shader.temperature = dummy_medium_getter;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &dummy));
+ return dummy;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)data; /* Avoid the "unused variable" warning */
+ return trilinear_profile(frag->P);
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+
+ /* Fluid/solid interface: fix the temperature to the temperature profile */
+ if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) {
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ }
+
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * Scene, i.e. the system to simulate
+ ******************************************************************************/
+struct scene_context {
+ const struct mesh* mesh;
+ size_t sshape_end_id; /* Last triangle index of the super shape */
+ struct sdis_interface* sshape;
+ struct sdis_interface* sphere;
+};
+#define SCENE_CONTEXT_NULL__ {NULL, 0, 0, 0}
+static const struct scene_context SCENE_CONTEXT_NULL = SCENE_CONTEXT_NULL__;
+
+static void
+scene_get_indices(const size_t itri, size_t ids[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(ids && context && itri < mesh_ntriangles(context->mesh));
+ /* Flip the indices to ensure that the normal points into the super shape */
+ ids[0] = (unsigned)context->mesh->indices[itri*3+0];
+ ids[1] = (unsigned)context->mesh->indices[itri*3+2];
+ ids[2] = (unsigned)context->mesh->indices[itri*3+1];
+}
+
+static void
+scene_get_interface(const size_t itri, struct sdis_interface** interf, void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(interf && context && itri < mesh_ntriangles(context->mesh));
+ if(itri < context->sshape_end_id) {
+ *interf = context->sshape;
+ } else {
+ *interf = context->sphere;
+ }
+}
+
+static void
+scene_get_position(const size_t ivert, double pos[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(pos && context && ivert < mesh_nvertices(context->mesh));
+ pos[0] = context->mesh->positions[ivert*3+0];
+ pos[1] = context->mesh->positions[ivert*3+1];
+ pos[2] = context->mesh->positions[ivert*3+2];
+}
+
+static struct sdis_scene*
+create_scene(struct sdis_device* sdis, struct scene_context* ctx)
+
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+
+ scn_args.get_indices = scene_get_indices;
+ scn_args.get_interface = scene_get_interface;
+ scn_args.get_position = scene_get_position;
+ scn_args.nprimitives = mesh_ntriangles(ctx->mesh);
+ scn_args.nvertices = mesh_nvertices(ctx->mesh);
+ scn_args.context = ctx;
+ OK(sdis_scene_create(sdis, &scn_args, &scn));
+ return scn;
+}
+
+/*******************************************************************************
+ * Validations
+ ******************************************************************************/
+static void
+check_probe_boundary_list_api
+ (struct sdis_scene* scn,
+ const struct sdis_solve_probe_boundary_list_args* in_args)
+{
+ struct sdis_solve_probe_boundary_list_args args = *in_args;
+ struct sdis_estimator_buffer* estim_buf = NULL;
+
+ /* Check API */
+ BA(sdis_solve_probe_boundary_list(NULL, &args, &estim_buf));
+ BA(sdis_solve_probe_boundary_list(scn, NULL, &estim_buf));
+ BA(sdis_solve_probe_boundary_list(scn, &args, NULL));
+ args.nprobes = 0;
+ BA(sdis_solve_probe_boundary_list(scn, &args, &estim_buf));
+ args.nprobes = in_args->nprobes;
+ args.probes = NULL;
+ BA(sdis_solve_probe_boundary_list(scn, &args, &estim_buf));
+}
+
+/* Check the estimators against the analytical solution */
+static void
+check_estimator_buffer
+ (const struct sdis_estimator_buffer* estim_buf,
+ struct sdis_scene* scn,
+ const struct sdis_solve_probe_boundary_list_args* args)
+{
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t iprobe = 0;
+
+ /* Variables used to check the global estimation results */
+ size_t total_nrealisations = 0;
+ size_t total_nrealisations_ref = 0;
+ size_t total_nfailures = 0;
+ size_t total_nfailures_ref = 0;
+ double sum = 0; /* Global sum of weights */
+ double sum2 = 0; /* Global sum of squared weights */
+ double N = 0; /* Number of (successful) realisations */
+ double E = 0; /* Expected value */
+ double V = 0; /* Variance */
+ double SE = 0; /* Standard Error */
+
+ /* Check the results */
+ FOR_EACH(iprobe, 0, args->nprobes) {
+ const struct sdis_estimator* estimator = NULL;
+ size_t probe_nrealisations = 0;
+ size_t probe_nfailures = 0;
+ double pos[3] = {0,0,0};
+ double probe_sum = 0;
+ double probe_sum2 = 0;
+ double ref = 0;
+
+ OK(sdis_scene_get_boundary_position(scn, args->probes[iprobe].iprim,
+ args->probes[iprobe].uv, pos));
+
+ /* Fetch result */
+ OK(sdis_estimator_buffer_at(estim_buf, iprobe, 0, &estimator));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_count(estimator, &probe_nrealisations));
+ OK(sdis_estimator_get_failure_count(estimator, &probe_nfailures));
+
+ /* Check probe estimation */
+ ref = trilinear_profile(pos);
+ printf("T(%g, %g, %g) = %g ~ %g +/- %g\n", SPLIT3(pos), ref, T.E, T.SE);
+ CHK(eq_eps(ref, T.E, 3*T.SE));
+
+ /* Check miscellaneous results */
+ CHK(probe_nrealisations == args->probes[iprobe].nrealisations);
+
+ /* Compute the sum of weights and sum of squared weights of the probe */
+ probe_sum = T.E * (double)probe_nrealisations;
+ probe_sum2 = (T.V + T.E*T.E) * (double)probe_nrealisations;
+
+ /* Update the global estimation */
+ total_nrealisations_ref += probe_nrealisations;
+ total_nfailures_ref += probe_nfailures;
+ sum += probe_sum;
+ sum2 += probe_sum2;
+ }
+
+ /* Check the overall estimate of the estimate buffer. This estimate is not
+ * really a result expected by the caller since the probes are all
+ * independent. But to be consistent with the estimate buffer API, we need to
+ * provide it. And so we check its validity */
+ OK(sdis_estimator_buffer_get_temperature(estim_buf, &T));
+ OK(sdis_estimator_buffer_get_realisation_count(estim_buf, &total_nrealisations));
+ OK(sdis_estimator_buffer_get_failure_count(estim_buf, &total_nfailures));
+
+ CHK(total_nrealisations == total_nrealisations_ref);
+ CHK(total_nfailures == total_nfailures_ref);
+
+ N = (double)total_nrealisations_ref;
+ E = sum / N;
+ V = sum2 / N - E*E;
+ SE = sqrt(V/N);
+ check_intersection(E, SE, T.E, T.SE);
+}
+
+static void
+check_probe_boundary_list(struct sdis_scene* scn, const int is_master_process)
+{
+ #define NPROBES 10
+
+ /* Estimations */
+ struct sdis_estimator_buffer* estim_buf = NULL;
+
+ /* Probe variables */
+ struct sdis_solve_probe_boundary_args probes[NPROBES];
+ struct sdis_solve_probe_boundary_list_args args =
+ SDIS_SOLVE_PROBE_BOUNDARY_LIST_ARGS_DEFAULT;
+ size_t iprobe;
+
+ (void)is_master_process;
+
+ /* Setup the list of probes to calculate */
+ args.probes = probes;
+ args.nprobes = NPROBES;
+ FOR_EACH(iprobe, 0, NPROBES) {
+ double pos[3] = {0, 0, 0};
+ sample_sphere(pos);
+
+ probes[iprobe] = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT;
+ probes[iprobe].nrealisations = 10000;
+ probes[iprobe].side = SDIS_FRONT;
+ OK(sdis_scene_find_closest_point
+ (scn, pos, INF, &probes[iprobe].iprim, probes[iprobe].uv));
+ }
+
+ check_probe_boundary_list_api(scn, &args);
+
+ /* Solve the probes */
+ OK(sdis_solve_probe_boundary_list(scn, &args, &estim_buf));
+ if(!is_master_process) {
+ CHK(estim_buf == NULL);
+ } else {
+ check_estimator_buffer(estim_buf, scn, &args);
+ OK(sdis_estimator_buffer_ref_put(estim_buf));
+ }
+
+ #undef NPROBES
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ /* Stardis */
+ struct sdis_device* dev = NULL;
+ struct sdis_interface* solid_fluid = NULL;
+ struct sdis_interface* solid_solid = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_scene* scn = NULL;
+
+ /* Miscellaneous */
+ struct scene_context ctx = SCENE_CONTEXT_NULL;
+ struct mesh mesh = MESH_NULL;
+ size_t sshape_end_id = 0; /* Last index of the super shape */
+ int is_master_process = 1;
+ (void)argc, (void)argv;
+
+ create_default_device(&argc, &argv, &is_master_process, &dev);
+
+ /* Setup the mesh */
+ mesh_init(&mesh);
+ mesh_add_super_shape(&mesh);
+ sshape_end_id = mesh_ntriangles(&mesh);
+ mesh_add_sphere(&mesh);
+
+ /* Setup physical properties */
+ fluid = create_dummy(dev);
+ solid = create_solid(dev);
+ solid_fluid = create_interface(dev, solid, fluid);
+ solid_solid = create_interface(dev, solid, solid);
+
+ /* Create the scene */
+ ctx.mesh = &mesh;
+ ctx.sshape_end_id = sshape_end_id;
+ ctx.sshape = solid_fluid;
+ ctx.sphere = solid_solid;
+ scn = create_scene(dev, &ctx);
+
+ check_probe_boundary_list(scn, is_master_process);
+
+ mesh_release(&mesh);
+ OK(sdis_interface_ref_put(solid_fluid));
+ OK(sdis_interface_ref_put(solid_solid));
+ OK(sdis_medium_ref_put(fluid));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_scene_ref_put(scn));
+
+ free_default_device(dev);
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -498,4 +498,3 @@ check_green_serialization
OK(sdis_estimator_ref_put(e2));
OK(sdis_green_function_ref_put(green2));
}
-
