commit 5463146a7990feda607c7388ffc795b5c06fa732
parent 8fe6ed7249c6ae9691099de1f2efb753867d2847
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 19 Mar 2026 11:57:25 +0100
sln-slab: improved error handling
Check that the importance sampling criterion is met at each stage of the
tree traversal, i.e., that the value of the node is greater than or equal
to the sum of the values of its children. If this is not the case,
display an error indicating issues related to the node values so that
the caller can verify the extent of the error and judge whether it is
more of an algorithm 'bug' or a numerical precision problem. The wave
number at which the error occurs is also displayed (with a sufficient
number of digits to ensure good precision), along with the path to the
node from the root, printed as a list of traversal options as defined by
the sln-get tool. This make easier the use of this tool to query this
node (and/or its children) at the wave number causing the problem, in
order to determine what is not working. This aims to improve the
debugging of the acceleration structure used to sample transitions.
As soon as an error is detected, the run fails, but the calculation
continues. The total number of rejected realisations is finally
displayed as a new output of the calculation.
Diffstat:
| M | src/sln_slab.c | | | 182 | ++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------- |
1 file changed, 121 insertions(+), 61 deletions(-)
diff --git a/src/sln_slab.c b/src/sln_slab.c
@@ -224,9 +224,11 @@ exit:
if(proxy) SSP(rng_proxy_ref_put(proxy));
return res;
error:
- fprintf(stderr,
- "Error creating the list of per thread RNG -- %s\n",
- res_to_cstr(res));
+ if(cmd->args.verbose >= 1) {
+ fprintf(stderr,
+ "Error creating the list of per thread RNG -- %s\n",
+ res_to_cstr(res));
+ }
if(rngs) delete_per_thread_rngs(cmd, rngs);
rngs = NULL;
goto exit;
@@ -297,13 +299,15 @@ error:
goto exit;
}
-static const struct sln_node*
+static const struct sln_node* /* NULL <=> an error occurs */
sample_line
(const struct cmd* cmd,
struct ssp_rng* rng,
const double nu/*[cm^-1]*/,
double* out_proba)
{
+ char step[64] = {0}; /* For debug */
+
const struct sln_node* node = NULL;
double proba = 1; /* Proba to sample the line */
size_t depth = 0;
@@ -314,9 +318,11 @@ sample_line
for(depth=0; ;++depth) {
const struct sln_node* node0 = NULL;
const struct sln_node* node1 = NULL;
+ struct sln_mesh mesh = SLN_MESH_NULL;
struct sln_mesh mesh0 = SLN_MESH_NULL;
struct sln_mesh mesh1 = SLN_MESH_NULL;
double proba0 = 0;
+ double ka_max = 0;
double ka_max0 = 0;
double ka_max1 = 0;
double r = 0;
@@ -331,15 +337,21 @@ sample_line
ka_max0 = sln_mesh_eval(&mesh0, nu);
ka_max1 = sln_mesh_eval(&mesh1, nu);
/* TODO handle ka_max[0,1] == 0 */
-#ifndef NDEBUG
- {
- struct sln_mesh mesh = SLN_MESH_NULL;
- double ka_max = 0;
- SLN(node_get_mesh(cmd->tree, node, &mesh));
- ka_max = sln_mesh_eval(&mesh, nu);
- ASSERT(ka_max >= ka_max0 + ka_max1);
+
+ /* Check the criterion of transition importance sampling, i.e. the value of
+ * the parent node is greater than or equal to the sum of the values of its
+ * children */
+ SLN(node_get_mesh(cmd->tree, node, &mesh));
+ ka_max = sln_mesh_eval(&mesh, nu);
+ if(ka_max < ka_max0 + ka_max1) {
+ step[depth] = '\0';
+ if(cmd->args.verbose >= 1) {
+ fprintf(stderr,
+ "error: ka < ka0 + ka1; %e < %e+%e; nu=%-21.20g cm^-1; node path: -%s\n",
+ ka_max, ka_max0, ka_max1, nu, step);
+ }
+ return NULL;
}
-#endif
/* Compute the probability to sample the 1st node */
proba0 = ka_max0 / (ka_max0 + ka_max1);
@@ -349,9 +361,11 @@ sample_line
if(r < proba0) {
node = node0;
proba *= proba0;
+ step[depth] = 'l';
} else {
node = node1;
proba *= (1-proba0);
+ step[depth] = 'r';
}
}
@@ -359,18 +373,58 @@ sample_line
return node;
}
-static double
-realisation(const struct cmd* cmd, struct ssp_rng* rng)
+static INLINE res_T
+check_sampled_node(const struct cmd* cmd, const struct sln_node* node)
+{
+ ASSERT(cmd);
+ (void) cmd;
+
+ if(!node) return RES_BAD_ARG;
+
+#ifndef NDEBUG
+ {
+ /* Verify that the sampled node corresponds to a single line */
+ struct sln_node_desc desc = SLN_NODE_DESC_NULL;
+ SLN(node_get_desc(node, &desc));
+ ASSERT(desc.nlines == 1);
+ }
+#endif
+
+ return RES_OK;
+}
+
+/* Check that the probability is valid */
+static INLINE res_T
+check_proba(const struct cmd* cmd, const double proba)
+{
+ if(proba <= 1) return RES_OK;
+
+ if(cmd->args.verbose >= 1) {
+ fprintf(stderr, "error: invalid probability %g\n", proba);
+ }
+ return RES_BAD_ARG;
+}
+
+static res_T
+realisation(const struct cmd* cmd, struct ssp_rng* rng, double* out_weight)
{
+ /* Acceleration structure */
const struct sln_node* node = NULL;
struct sln_mesh mesh = SLN_MESH_NULL;
+
+ /* Null collisions */
double ka_max = 0;
double dst_remain = 0;
- double dst = 0; /* Length */
- double nu = 0;
- size_t ncollisions = 0; /* For debug */
- ASSERT(cmd && rng);
+ size_t ncollisions = 0; /* Number of null collisions */
+
+ /* Miscellaneous */
+ double nu = 0; /* Sampled wavenumber */
+ double w = 0; /* Monte Carlo weight */
+ res_T res = RES_OK;
+
+ ASSERT(cmd && rng && out_weight); /* Check pre-conditions */
+ /* Initialize the total distance to traverse with the thickness of the slab */
dst_remain = cmd->args.thickness;
/* Uniformly sample the spectral dimension */
@@ -383,45 +437,41 @@ realisation(const struct cmd* cmd, struct ssp_rng* rng)
ka_max = sln_mesh_eval(&mesh, nu);
for(ncollisions=0; ; ++ncollisions) {
- /* Sample a traversal distance */
- dst = ssp_ran_exp(rng, ka_max);
+ double dst = 0; /* Sampled distance */
+ double proba_abs = 0; /* Probability of absorption */
+ double line_proba = 0; /* Probability of sampling the line */
+ double line_ha = 0; /* Value of the line */
+ double r = 0; /* Random number */
- if(dst > dst_remain) {
- return cmd->args.spectral_range[1] - cmd->args.spectral_range[0];
- } else {
- double proba_abs = 0;
- double line_proba = 0;
- double line_ha = 0; /* Value of the line */
- double r = 0; /* Random number */
-
- /* Importance sampling of a line */
- node = sample_line(cmd, rng, nu, &line_proba);
-
- /* Check whether the sampled node corresponds to a sampled line */
- #ifndef NDEBUG
- {
- struct sln_node_desc desc = SLN_NODE_DESC_NULL;
- SLN(node_get_desc(node, &desc));
- ASSERT(desc.nlines == 1);
- }
- #endif
+ dst = ssp_ran_exp(rng, ka_max); /* Sample a traversal distance */
- /* Get the line value */
- line_ha = sln_node_eval(cmd->tree, node, nu);
+ if(dst > dst_remain) { /* No absorption in the slab */
+ w = cmd->args.spectral_range[1] - cmd->args.spectral_range[0];
+ break;
+ }
- proba_abs = line_ha / (line_proba*ka_max);
- if(proba_abs > 1) {
- VFATAL("Invalid probability %g\n", ARG1(proba_abs));
- }
+ /* Importance sampling of a line */
+ node = sample_line(cmd, rng, nu, &line_proba);
+ if((res = check_sampled_node(cmd, node)) != RES_OK) goto error;
- r = ssp_rng_canonical(rng);
- if(r < proba_abs) {
- return 0;
- } else { /* Null collision */
- dst_remain -= dst;
- }
- }
+ /* Evaluate the value of the line and compute the probability of being
+ * absorbed by it */
+ line_ha = sln_node_eval(cmd->tree, node, nu);
+ proba_abs = line_ha / (line_proba*ka_max);
+ if((res = check_proba(cmd, proba_abs)) != RES_OK) goto error;
+
+ r = ssp_rng_canonical(rng);
+ if(r < proba_abs) { w=0; break; } /* An absorption occurs */
+
+ dst_remain -= dst; /* This was a null transition. Go on */
}
+
+exit:
+ *out_weight = w;
+ return res;
+error:
+ w = NaN;
+ goto exit;
}
static res_T
@@ -436,6 +486,7 @@ cmd_run(const struct cmd* cmd)
double V = 0; /* Variance */
double SE = 0; /* Standard Error */
int64_t i = 0; /* Index of the realisation */
+ size_t nrejects = 0; /* Number of rejected realisations */
/* Progress */
size_t nrealisations = 0;
@@ -449,42 +500,51 @@ cmd_run(const struct cmd* cmd)
res = create_per_thread_rngs(cmd, &rngs);
if(res != RES_OK) goto error;
- if(cmd->args.verbose) fprintf(stderr, "%3d%%\n", progress);
+ #define PROGRESS_MSG "Solving: %3d%%\n"
+ if(cmd->args.verbose >= 3) fprintf(stderr, PROGRESS_MSG, progress);
nrealisations = cmd->args.nrealisations;
+ omp_set_num_threads((int)cmd->nthreads);
+
#pragma omp parallel for schedule(static)
for(i = 0; i < (int64_t)nrealisations; ++i) {
const int ithread = omp_get_thread_num();
int pcent = 0;
- double w = 0;
+ double w = 0; /* Monte Carlo weight */
+ res_T res_realisation = RES_OK;
- w = realisation(cmd, rngs[ithread]);
+ res_realisation = realisation(cmd, rngs[ithread], &w);
#pragma omp critical
{
/* Update the Monte Carlo accumulator */
- accum.sum += w;
- accum.sum2 += w*w;
- accum.count += 1;
+ if(res_realisation == RES_OK) {
+ accum.sum += w;
+ accum.sum2 += w*w;
+ accum.count += 1;
+ }
- if(cmd->args.verbose) {
+ if(cmd->args.verbose >= 3) {
/* Update progress */
realisation_done += 1;
pcent = (int)((double)realisation_done*100.0/(double)nrealisations+0.5);
if(pcent/progress_pcent > progress/progress_pcent) {
progress = pcent;
- fprintf(stderr, "%3d%%\n", progress);
+ fprintf(stderr, PROGRESS_MSG, progress);
}
}
}
}
+ #undef PROGRESS_MSG
+
E = accum.sum / (double)accum.count;
V = accum.sum2 / (double)accum.count - E*E;
SE = sqrt(V/(double)accum.count);
+ nrejects = nrealisations - accum.count;
- printf("%g %g\n", E, SE);
+ printf("%g %g %lu\n", E, SE, (unsigned long)nrejects);
exit:
delete_per_thread_rngs(cmd, rngs);
