commit 90a512e2d9cd6b7b56f03574c40cfba31db4dee6
parent 9468a804126acd6ba536ba8894c5532ef957bace
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 12 Mar 2026 16:33:04 +0100
sln-slab: preparation for Monte Carlo simulation
Setup the parallel estimation of the Monte Carlo algorithm, which
currently calculates 4+1. With a verbosity of 3, a computation progress
message is displayed.
The next step is to implement the actual function for performing the
radiative transfer algorithm.
Diffstat:
| M | Makefile | | | 30 | ++++++++++++++++++++++++------ |
| M | config.mk | | | 1 | + |
| M | src/sln_slab.c | | | 142 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
3 files changed, 167 insertions(+), 6 deletions(-)
diff --git a/Makefile b/Makefile
@@ -70,6 +70,7 @@ libsln.o: $(OBJ)
$(PKG_CONFIG) --atleast-version $(LBLU_VERSION) lblu
$(PKG_CONFIG) --atleast-version $(RSYS_VERSION) rsys
$(PKG_CONFIG) --atleast-version $(SHTR_VERSION) shtr
+ $(PKG_CONFIG) --atleast-version $(SSP_VERSION) star-sp
echo "config done" > $@
.SUFFIXES: .c .d .o
@@ -90,15 +91,26 @@ PKG_CONFIG_LOCAL = PKG_CONFIG_PATH="./:$${PKG_CONFIG_PATH}" $(PKG_CONFIG)
INCS_UTIL = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --cflags rsys sln-local shtr)
LIBS_UTIL = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys sln-local shtr)
-
CFLAGS_UTIL = -std=c89 $(CFLAGS_EXE) $(INCS_UTIL)
LDFLAGS_UTIL = $(LDFLAGS_EXE) $(LIBS_UTIL)
-utils: library $(UTIL_DEP)
+INCS_SLAB = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --cflags rsys sln-local shtr star-sp)
+LIBS_SLAB = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys sln-local shtr star-sp) -lm
+CFLAGS_SLAB = -std=c89 $(CFLAGS_EXE) $(INCS_SLAB) -fopenmp
+LDFLAGS_SLAB = $(LDFLAGS_EXE) $(LIBS_SLAB)
+
+utils: library $(UTIL_DEP) .config
@$(MAKE) -fMakefile \
$$(for i in $(UTIL_DEP); do printf -- '-f%s\n' "$${i}"; done) \
sln-build sln-get sln-slab
+$(UTIL_DEP) $(UTIL_OBJ): config.mk sln-local.pc
+
+
+src/sln_build.d src/sln_build.o: src/sln_build.c
+src/sln_get.d src/sln_get.o: src/sln_get.c
+src/sln_slab.d src/sln_slab.o: src/sln_slab.c
+
sln-build: config.mk sln-local.pc src/sln_build.o $(LIBNAME)
$(CC) $(CFLAGS_UTIL) -o $@ src/sln_build.o $(LDFLAGS_UTIL)
@@ -106,16 +118,22 @@ sln-get: config.mk sln-local.pc src/sln_get.o $(LIBNAME)
$(CC) $(CFLAGS_UTIL) -o $@ src/sln_get.o $(LDFLAGS_UTIL)
sln-slab: config.mk sln-local.pc src/sln_slab.o $(LIBNAME)
- $(CC) $(CFLAGS_UTIL) -o $@ src/sln_slab.o $(LDFLAGS_UTIL)
+ $(CC) $(CFLAGS_SLAB) -o $@ src/sln_slab.o $(LDFLAGS_SLAB)
-$(UTIL_DEP): config.mk sln-local.pc
+src/sln_build.d src/sln_get.d:
@$(CC) $(CFLAGS_UTIL) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
-$(UTIL_OBJ): config.mk sln-local.pc
+src/sln_slab.d:
+ @$(CC) $(CFLAGS_SLAB) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
+
+src/sln_build.o src/sln_get.o:
$(CC) $(CFLAGS_UTIL) -c $(@:.o=.c) -o $@
+src/sln_slab.o:
+ $(CC) $(CFLAGS_SLAB) -c $(@:.o=.c) -o $@
+
clean_utils:
- rm -f $(UTIL_OBJ) $(UTIL_DEP) sln-build
+ rm -f $(UTIL_OBJ) $(UTIL_DEP) sln-build sln-get sln-slab
################################################################################
# Installation
diff --git a/config.mk b/config.mk
@@ -31,6 +31,7 @@ PCFLAGS = $(PCFLAGS_$(LIB_TYPE))
LBLU_VERSION = 0.0.1
RSYS_VERSION = 0.14
SHTR_VERSION = 0
+SSP_VERSION = 0.15 # Dependency for the sln-slab util
INCS = $$($(PKG_CONFIG) $(PCFLAGS) --cflags lblu rsys shtr)
LIBS = $$($(PKG_CONFIG) $(PCFLAGS) --libs lblu rsys shtr) -lm
diff --git a/src/sln_slab.c b/src/sln_slab.c
@@ -21,10 +21,13 @@
#include "sln.h"
#include <star/shtr.h>
+#include <star/ssp.h>
#include <rsys/cstr.h>
#include <rsys/mem_allocator.h>
+#include <omp.h>
+
#include <unistd.h> /* getopt */
struct args {
@@ -49,10 +52,19 @@ struct cmd {
struct args args;
struct sln_tree* tree;
+ unsigned nthreads;
};
#define CMD_NULL__ {0}
static const struct cmd CMD_NULL = CMD_NULL__;
+struct accum {
+ double sum;
+ double sum2;
+ size_t count;
+};
+#define ACCUM_NULL__ {0}
+static const struct accum ACCUM_NULL = ACCUM_NULL__;
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -159,6 +171,50 @@ error:
goto exit;
}
+static void
+delete_per_thread_rngs(const struct cmd* cmd, struct ssp_rng* rngs[])
+{
+ unsigned i = 0;
+ ASSERT(cmd && rngs);
+
+ FOR_EACH(i, 0, cmd->nthreads) {
+ if(rngs[i]) SSP(rng_ref_put(rngs[i]));
+ }
+ mem_rm(rngs);
+}
+
+static res_T
+create_per_thread_rngs(const struct cmd* cmd, struct ssp_rng** out_rngs[])
+{
+ struct ssp_rng_proxy* proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ size_t i = 0;
+ res_T res = RES_OK;
+ ASSERT(cmd);
+
+ rngs = mem_calloc(cmd->nthreads, sizeof(*rngs));
+ if(!rngs) { res = RES_MEM_ERR; goto error; }
+
+ res = ssp_rng_proxy_create(NULL, SSP_RNG_THREEFRY, cmd->nthreads, &proxy);
+ if(res != RES_OK) goto error;
+
+ FOR_EACH(i, 0, cmd->nthreads) {
+ res = ssp_rng_proxy_create_rng(proxy, i, &rngs[i]);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ *out_rngs = rngs;
+ 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(rngs) delete_per_thread_rngs(cmd, rngs);
+ rngs = NULL;
+ goto exit;
+}
static void
cmd_release(struct cmd* cmd)
@@ -182,6 +238,7 @@ cmd_init(struct cmd* cmd, const struct args* args)
struct shtr_line_list* lines = NULL;
/* Miscellaneous */
+ unsigned nthreads_max = 0;
res_T res = RES_OK;
ASSERT(cmd && args);
@@ -208,7 +265,9 @@ cmd_init(struct cmd* cmd, const struct args* args)
res = sln_tree_read(sln, &tree_args, &cmd->tree);
if(res != RES_OK) goto error;
+ nthreads_max = (unsigned)MMAX(omp_get_max_threads(), omp_get_num_procs());
cmd->args = *args;
+ cmd->nthreads = MMIN(cmd->args.nthreads_hint, nthreads_max);
exit:
if(sln) SLN(device_ref_put(sln));
@@ -222,6 +281,88 @@ error:
goto exit;
}
+static double
+realisation(const struct cmd* cmd, struct ssp_rng* rng)
+{
+ ASSERT(cmd && rng);
+ (void)cmd; /* Avoid unused variable warning */
+
+ if(ssp_rng_canonical(rng) < 0.5) {
+ return 2.0;
+ } else {
+ return 8.0;
+ }
+}
+
+static res_T
+cmd_run(const struct cmd* cmd)
+{
+ /* Random Number Generator */
+ struct ssp_rng** rngs = NULL;
+
+ /* Monte Carlo */
+ struct accum accum = ACCUM_NULL;
+ double E = 0; /* Expected value */
+ double V = 0; /* Variance */
+ double SE = 0; /* Standard Error */
+ int64_t i = 0; /* Index of the realisation */
+
+ /* Progress */
+ size_t nrealisations = 0;
+ size_t realisation_done = 0;
+ int progress = 0;
+ int progress_pcent = 10;
+
+ res_T res = RES_OK;
+ ASSERT(cmd);
+
+ res = create_per_thread_rngs(cmd, &rngs);
+ if(res != RES_OK) goto error;
+
+ if(cmd->args.verbose) fprintf(stderr, "%3d%%\n", progress);
+
+ nrealisations = cmd->args.nrealisations;
+
+ #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;
+
+ w = realisation(cmd, rngs[ithread]);
+
+ #pragma omp critical
+ {
+ /* Update the Monte Carlo accumulator */
+ accum.sum += w;
+ accum.sum2 += w*w;
+ accum.count += 1;
+
+ if(cmd->args.verbose) {
+ /* 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);
+ }
+ }
+ }
+ }
+
+ E = accum.sum / (double)accum.count;
+ V = accum.sum2 / (double)accum.count - E*E;
+ SE = sqrt(V/(double)accum.count);
+
+ printf("%g %g\n", E, SE);
+
+exit:
+ delete_per_thread_rngs(cmd, rngs);
+ return res;
+error:
+ goto exit;
+}
+
/*******************************************************************************
* Main function
******************************************************************************/
@@ -237,6 +378,7 @@ main(int argc, char** argv)
if(args.quit) goto exit;
if((res = cmd_init(&cmd, &args)) != RES_OK) goto error;
+ if((res = cmd_run(&cmd)) != RES_OK) goto error;
exit:
cmd_release(&cmd);
