32  Minimal HPC simulation template (Everything)

33 What this is

A stripped-down template for running a parallel Monte Carlo simulation on the “Everything” HPC (SLURM + interactive RStudio Server, ~196 cores but 4-16 is the max you should use; little performance gains beyond this) that also runs locally without changes. It keeps only the cluster scaffolding — worker auto-detection, BLAS thread pinning, reproducible parallel seeding, condition-level parallelism with summarise-in-worker, and caching — wrapped around a single trivial simulation: the statistical power of a two-sample t-test. Swap the simulate_trial() and run_condition() functions for your own; the orchestration below is the reusable part.

The design choices that matter on a cluster:

• Parallelise over conditions, not replicates. Each worker runs all replicates for one condition and reduces them to a few summary numbers before returning, so the master process never holds the raw replicate rows. Memory stays flat regardless of how large the factorial gets.
• One replication knob, auto-scaled. The same file does a fast local smoke test and the full cluster run; on_hpc detects which and sets n_iterations. Caches are tagged so the two never collide.
• Cache the summaries, not the raw draws. Re-rendering the report reads a small .rds instead of recomputing.

34 Dependencies and parallel plan

library(tidyr)
library(dplyr)
library(furrr)
library(readr)
library(ggplot2)

# --- How many workers does THIS session/allocation actually have? -----------
# parallelly::availableCores() honours the SLURM allocation, cgroup CPU quotas,
# and affinity masks, so it returns the cores granted to THIS session -- not the
# ~196 physical cores parallel::detectCores() would report. Set the override to
# pin a count manually; leave NULL to auto-detect.
n_workers_override <- NULL
detected_cores <- unname(parallelly::availableCores())
n_workers <- if (is.null(n_workers_override)) detected_cores else as.integer(n_workers_override)
n_workers <- max(1L, as.integer(n_workers))
message(sprintf("Parallel plan: %d workers (parallelly detected %d).", n_workers, detected_cores))

# --- Pin BLAS/LAPACK to one thread BEFORE workers spawn. --------------------
# Everything's R uses FlexiBLAS over OpenBLAS, which defaults to many threads
# (blas_get_num_procs() reported 128). With dozens of worker processes that is
# 64 x 128 threads fighting over the node. escalc()/rma() use tiny matrices, so
# single-threaded BLAS loses nothing. Setting these in the master means the
# multisession (PSOCK) workers inherit them at spawn.
Sys.setenv(FLEXIBLAS_NUM_THREADS = "1", OMP_NUM_THREADS = "1",
           OPENBLAS_NUM_THREADS = "1", MKL_NUM_THREADS = "1", BLIS_NUM_THREADS = "1",
           VECLIB_MAXIMUM_THREADS = "1")  # VECLIB covers macOS Accelerate (local testing)
if (requireNamespace("flexiblas", quietly = TRUE)) flexiblas::flexiblas_set_num_threads(1L)
if (requireNamespace("RhpcBLASctl", quietly = TRUE)) {
  RhpcBLASctl::blas_set_num_threads(1)
  RhpcBLASctl::omp_set_num_threads(1)
}

# RStudio Server is not fork-safe, so use multisession (separate R processes),
# not multicore. furrr_options(seed = TRUE) (below) advances an L'Ecuyer-CMRG
# stream per task, giving reproducibility independent of worker count.
plan(multisession, workers = n_workers)

35 Replication knob (local vs cluster)

# Detection note: RStudio Server does NOT pass the job's SLURM_* allocation
# variables into the R session, but Everything sets SLURM_MPI_TYPE as a global
# default that DOES survive and is absent on a laptop -- so "any SLURM_* variable
# present" is a reliable cluster detector (more robust than core count, since an
# allocation may be as small as 4 cores). run_full forces the choice manually.
run_full <- NULL  # NULL = auto-detect; TRUE = full run; FALSE = quick local test

on_hpc <- if (!is.null(run_full)) isTRUE(run_full) else any(grepl("^SLURM_", names(Sys.getenv())))
n_iterations <- if (on_hpc) 1000L else 50L

36 The simulation (swap these for your own)

Two small functions: one replicate (simulate two groups and test them), and the condition-level reduction that runs many replicates and collapses them to a power estimate. Only the second is HPC-specific.

# --- One replicate: two independent groups of size n, separated by delta. ----
# Returns TRUE if the two-sample t-test rejects at the given alpha.
simulate_trial <- function(n, delta, alpha = 0.05) {
  x <- rnorm(n, mean = 0)
  y <- rnorm(n, mean = delta)
  t.test(x, y)$p.value < alpha
}

# --- Condition-level reduction (the HPC-relevant part). ----------------------
# Run all n_iter replicates for ONE condition, then collapse to a single
# performance number -- the rejection rate (power; or Type I error when
# delta = 0) -- with its Monte Carlo standard error, before returning. The
# worker hands back a couple of numbers, never the raw replicate results.
run_condition <- function(n, delta, alpha, n_iter) {
  rejects <- replicate(n_iter, simulate_trial(n, delta, alpha))
  power   <- mean(rejects)
  tibble(n_iter = n_iter,
         power = power,
         power_mcse = sqrt(power * (1 - power) / n_iter))   # SE of a proportion
}

37 Experiment grid

# Locally, collapse to a small slice so a render takes seconds; on the HPC use
# the full factorial.
.n <- if (on_hpc) c(10, 20, 30, 50, 75, 100) else c(20, 100)

# set seed for load balancing randomisation via slice_ below, to ensure subsequent results are reproducible
set.seed(42) 

experiment_parameters <- expand_grid(
  n     = .n,                              # per-group sample size
  delta = c(0, 0.2, 0.5, 0.8),             # true mean difference (delta = 0 -> Type I rate)
  alpha = 0.05
) |>
  # shuffle the rows of the grid. 
  # this often speeds up simulations, as it is a form of 'load balancing': when being executed in parallel, it ensures that the computationally expensive iterations (eg with highest N) are distributed across cores/workers
  slice_sample(prop = 1) 

38 Run (cached)

cache_file <- paste0("simulation_minimal_", if (on_hpc) "full" else "dev",
                     "_k", n_iterations, ".rds")

if (file.exists(cache_file)) {
  simulation <- read_rds(cache_file)
} else {
  set.seed(42)
  simulation <- experiment_parameters |>
    mutate(summary = future_pmap(
      .l = list(n = n, delta = delta, alpha = alpha),
      .f = run_condition,
      n_iter = n_iterations,
      .progress = TRUE,
      .options = furrr_options(
        seed = TRUE, scheduling = 2,
        packages = c("dplyr", "tibble"))
    )) |>
    unnest(summary) |>
    arrange(delta, n)
  write_rds(simulation, cache_file, compress = "gz")
}

39 A look at the result

simulation |>
  ggplot(aes(n, power, colour = factor(delta))) +
  geom_hline(yintercept = c(0.05, 0.80), linetype = "dashed", colour = "grey70") +
  geom_line() +
  geom_point(size = 1.8) +
  scale_colour_viridis_d(name = "True delta", end = 0.85) +
  scale_x_continuous(name = "Per-group sample size n") +
  scale_y_continuous(name = "Rejection rate (power)", limits = c(0, 1)) +
  theme_linedraw(base_size = 12) +
  theme(panel.grid.minor = element_blank(), legend.position = "top") +
  ggtitle("Power of a two-sample t-test",
          subtitle = "delta = 0 traces the Type I error rate (~0.05); dashed lines at 0.05 and 0.80")

40 Running on Everything (SLURM + RStudio Server)

• Request a node and launch RStudio Server inside the allocation, e.g. salloc --nodes=1 --cpus-per-task=64 --mem=128G --time=04:00:00. The config chunk auto-detects the granted cores via parallelly::availableCores() and adapts.
• Belt-and-braces for BLAS: export the thread limit in the shell before launching R — export FLEXIBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 — and confirm with RhpcBLASctl::blas_get_num_procs() (should print 1).
• multisession, not multicore: RStudio is not fork-safe. If you run headless via Rscript/quarto render on Linux, multicore workers are cheaper.
• Reproducible regardless of worker count: each condition gets its own L’Ecuyer stream by task index, so changing n_workers does not change the numbers.
• Single-node (shared-memory) parallelism. Spanning nodes would need a different backend (e.g. future.batchtools with a SLURM template) and is unnecessary for this workload.

sessionInfo()

R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.5

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Zurich
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ggplot2_4.0.3 readr_2.2.0   furrr_0.3.1   future_1.67.0 dplyr_1.2.1  
[6] tidyr_1.3.2  

loaded via a namespace (and not attached):
 [1] gtable_0.3.6        jsonlite_2.0.0      compiler_4.5.2     
 [4] tidyselect_1.2.1    parallel_4.5.2      scales_1.4.0       
 [7] globals_0.18.0      RhpcBLASctl_0.23-42 yaml_2.3.12        
[10] fastmap_1.2.0       R6_2.6.1            labeling_0.4.3     
[13] flexiblas_3.4.0     generics_0.1.4      knitr_1.51         
[16] htmlwidgets_1.6.4   tibble_3.3.1        RColorBrewer_1.1-3 
[19] pillar_1.11.1       tzdb_0.5.0          rlang_1.3.0        
[22] xfun_0.60           S7_0.2.2            otel_0.2.0         
[25] viridisLite_0.4.3   cli_3.6.6           withr_3.0.3        
[28] magrittr_2.0.5      digest_0.6.39       grid_4.5.2         
[31] rstudioapi_0.19.0   hms_1.1.4           lifecycle_1.0.5    
[34] vctrs_0.7.3         evaluate_1.0.5      glue_1.8.1         
[37] farver_2.1.2        listenv_0.9.1       codetools_0.2-20   
[40] parallelly_1.45.1   rmarkdown_2.31      purrr_1.2.2        
[43] tools_4.5.2         pkgconfig_2.0.3     htmltools_0.5.9