commit 2905793c247efa1ae8b7c8e60ff16fb568a71252
parent 2c4a5b8b76f79995692698f70eea9cf9aa6a224e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 11 Dec 2023 16:45:27 +0100
Further tests on the resolution of a list of probes
Test the global results returned by the estimator buffer and check that
the state of the input RNG is correctly taken into account.
Diffstat:
3 files changed, 234 insertions(+), 24 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -9,5 +9,6 @@ test*
.config
.config_test
.test
+rng_state
tags
*.pc
diff --git a/Makefile b/Makefile
@@ -255,6 +255,7 @@ test_all: test
clean_test:
@$(SHELL) make.sh clean_test $(TEST_SRC) $(TEST_SRC_MPI) $(TEST_SRC_LONG)
+ rm -f rng_state
################################################################################
# Regular tests
diff --git a/src/test_sdis_solve_probe_list.c b/src/test_sdis_solve_probe_list.c
@@ -20,6 +20,11 @@
#include <star/s3d.h>
#include <star/s3dut.h>
+#include <star/ssp.h>
+
+#ifndef SDIS_ENABLE_MPI
+ #include <mpi.h>
+#endif
/*
* The system is a trilinear profile of steady state temperature, i.e. at each
@@ -75,6 +80,24 @@ rand_canonic(void)
return (float)rand() / (float)((unsigned)RAND_MAX+1);
}
+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]);
+}
+
/*******************************************************************************
* Super shape
******************************************************************************/
@@ -372,6 +395,161 @@ create_interface
* Validations
******************************************************************************/
static void
+check_probe_list_api
+ (struct sdis_scene* scn,
+ const struct sdis_solve_probe_list_args* in_args)
+{
+ struct sdis_solve_probe_list_args args = *in_args;
+ struct sdis_estimator_buffer* estim_buf = NULL;
+
+ /* Check API */
+ BA(sdis_solve_probe_list(NULL, &args, &estim_buf));
+ BA(sdis_solve_probe_list(scn, NULL, &estim_buf));
+ BA(sdis_solve_probe_list(scn, &args, NULL));
+ args.nprobes = 0;
+ BA(sdis_solve_probe_list(scn, &args, &estim_buf));
+ args.nprobes = in_args->nprobes;
+ args.probes = NULL;
+ BA(sdis_solve_probe_list(scn, &args, &estim_buf));
+}
+
+/* Check the estimators against the analytical solution */
+static void
+check_estimator_buffer
+ (const struct sdis_estimator_buffer* estim_buf,
+ const struct sdis_solve_probe_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 probe_sum = 0;
+ double probe_sum2 = 0;
+ double ref = 0;
+
+ /* 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(args->probes[iprobe].position);
+ printf("T(%g, %g, %g) = %g ~ %g +/- %g\n",
+ SPLIT3(args->probes[iprobe].position), 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);
+}
+
+/* Check that the buffers store statistically independent estimates */
+static void
+check_estimator_buffer_combination
+ (const struct sdis_estimator_buffer* estim_buf1,
+ const struct sdis_estimator_buffer* estim_buf2,
+ const struct sdis_solve_probe_list_args* args)
+{
+ size_t iprobe;
+
+ /* Check that the 2 series of estimates are compatible but not identical */
+ FOR_EACH(iprobe, 0, args->nprobes) {
+ const struct sdis_estimator* estimator1 = NULL;
+ const struct sdis_estimator* estimator2 = NULL;
+ size_t nrealisations1 = 0;
+ size_t nrealisations2 = 0;
+ struct sdis_mc T1 = SDIS_MC_NULL;
+ struct sdis_mc T2 = SDIS_MC_NULL;
+ double sum = 0; /* Sum of weights */
+ double sum2 = 0; /* Sum of squared weights */
+ double E = 0; /* Expected value */
+ double V = 0; /* Variance */
+ double SE = 0; /* Standard Error */
+ double N = 0; /* Number of realisations */
+ double ref = 0; /* Analytical solution */
+
+ /* Retrieve the estimation results */
+ OK(sdis_estimator_buffer_at(estim_buf1, iprobe, 0, &estimator1));
+ OK(sdis_estimator_buffer_at(estim_buf2, iprobe, 0, &estimator2));
+ OK(sdis_estimator_get_temperature(estimator1, &T1));
+ OK(sdis_estimator_get_temperature(estimator2, &T2));
+
+ /* Estimates are not identical... */
+ CHK(T1.E != T2.E);
+ CHK(T1.SE != T2.SE);
+
+ /* ... but are compatible ... */
+ check_intersection(T1.E, 3*T1.SE, T2.E, 3*T2.SE);
+
+ /* We can combine the 2 estimates since they must be numerically
+ * independent. We can therefore verify that their combination is valid with
+ * respect to the analytical solution */
+ OK(sdis_estimator_get_realisation_count(estimator1, &nrealisations1));
+ OK(sdis_estimator_get_realisation_count(estimator2, &nrealisations2));
+
+ sum =
+ T1.E*(double)nrealisations1
+ + T2.E*(double)nrealisations2;
+ sum2 =
+ (T1.V + T1.E*T1.E)*(double)nrealisations1
+ + (T2.V + T2.E*T2.E)*(double)nrealisations2;
+ N = (double)(nrealisations1 + nrealisations2);
+ E = sum / N;
+ V = sum2 / N - E*E;
+ SE = sqrt(V/N);
+
+ ref = trilinear_profile(args->probes[iprobe].position);
+ printf("T(%g, %g, %g) = %g ~ %g +/- %g\n",
+ SPLIT3(args->probes[iprobe].position), ref, E, SE);
+ CHK(eq_eps(ref, E, 3*SE));
+ }
+}
+
+static void
check_probe_list
(struct sdis_scene* scn,
struct s3d_scene_view* view,
@@ -379,16 +557,23 @@ check_probe_list
{
#define NPROBES 10
+ /* RNG state */
+ const char* rng_state_filename = "rng_state";
+ struct ssp_rng* rng = NULL;
+ enum ssp_rng_type rng_type = SSP_RNG_TYPE_NULL;
+
/* Probe variables */
struct sdis_solve_probe_args probes[NPROBES];
struct sdis_solve_probe_list_args args = SDIS_SOLVE_PROBE_LIST_ARGS_DEFAULT;
size_t iprobe = 0;
- struct sdis_mc T = SDIS_MC_NULL;
+ /* Estimations */
struct sdis_estimator_buffer* estim_buf = NULL;
+ struct sdis_estimator_buffer* estim_buf2 = NULL;
- double ref = 0; /* Analytical reference */
+ /* Misc */
double delta = 0;
+ FILE* fp = NULL;
delta = view_compute_delta(view);
@@ -401,38 +586,61 @@ check_probe_list
view_sample_position(view, delta, probes[iprobe].position);
}
- /* Check API */
- BA(sdis_solve_probe_list(NULL, &args, &estim_buf));
- BA(sdis_solve_probe_list(scn, NULL, &estim_buf));
- BA(sdis_solve_probe_list(scn, &args, NULL));
- args.nprobes = 0;
- BA(sdis_solve_probe_list(scn, &args, &estim_buf));
- args.nprobes = NPROBES;
- args.probes = NULL;
- BA(sdis_solve_probe_list(scn, &args, &estim_buf));
- args.probes = probes;
+ check_probe_list_api(scn, &args);
/* Solve the probes */
OK(sdis_solve_probe_list(scn, &args, &estim_buf));
if(!is_master_process) {
CHK(estim_buf == NULL);
- return; /* Nothing to do */
+ } else {
+ check_estimator_buffer(estim_buf, &args);
+
+ /* Check the use of a non-default rng_state. Run the calculations using the
+ * final RNG state from the previous estimate. Thus, the estimates are
+ * statically independent and can be combined. To do this, write the RNG
+ * state to a file so that it is read by all processes for the next test. */
+ OK(sdis_estimator_buffer_get_rng_state(estim_buf, &rng));
+ OK(ssp_rng_get_type(rng, &rng_type));
+ CHK(fp = fopen(rng_state_filename, "w"));
+ CHK(fwrite(&rng_type, sizeof(rng_type), 1, fp) == 1);
+ OK(ssp_rng_write(rng, fp));
+ CHK(fclose(fp) == 0);
}
- /* Check the results */
- FOR_EACH(iprobe, 0, NPROBES) {
- const struct sdis_estimator* estimator = NULL;
- OK(sdis_estimator_buffer_at(estim_buf, iprobe, 0, &estimator));
- OK(sdis_estimator_get_temperature(estimator, &T));
-
- ref = trilinear_profile(probes[iprobe].position);
- printf("T(%g, %g, %g) = %g ~ %g +/- %g\n",
- SPLIT3(probes[iprobe].position), ref, T.E, T.SE);
- CHK(eq_eps(ref, T.E, 3*T.SE));
+#ifdef SDIS_ENABLE_MPI
+ /* Synchronize processes. Wait for the master process to write RNG state to
+ * be used */
+ {
+ struct sdis_device* sdis = NULL;
+ int is_mpi_used = 0;
+ OK(sdis_scene_get_device(scn, &sdis));
+ OK(sdis_device_is_mpi_used(sdis, &is_mpi_used));
+ if(is_mpi_used) {
+ CHK(MPI_Barrier(MPI_COMM_WORLD) == MPI_SUCCESS);
+ }
}
+#endif
- OK(sdis_estimator_buffer_ref_put(estim_buf));
+ /* Read the RNG state */
+ CHK(fp = fopen(rng_state_filename, "r"));
+ CHK(fread(&rng_type, sizeof(rng_type), 1, fp) == 1);
+ OK(ssp_rng_create(NULL, rng_type, &rng));
+ OK(ssp_rng_read(rng, fp));
+ CHK(fclose(fp) == 0);
+
+ /* Run a new calculation using the read RNG state */
+ args.rng_state = rng;
+ OK(sdis_solve_probe_list(scn, &args, &estim_buf2));
+ OK(ssp_rng_ref_put(rng));
+
+ if(!is_master_process) {
+ CHK(estim_buf2 == NULL);
+ } else {
+ check_estimator_buffer_combination(estim_buf, estim_buf2, &args);
+ OK(sdis_estimator_buffer_ref_put(estim_buf));
+ OK(sdis_estimator_buffer_ref_put(estim_buf2));
+ }
#undef NPROBES
}
