From f8ecbae99dd11097c18ebfec361d6b1632415756 Mon Sep 17 00:00:00 2001
From: salman <zaher.salman@psi.ch>
Date: Sun, 14 Jun 2026 22:28:29 +0200
Subject: [PATCH] Make scans run on multiple threads. Speedup gain on local
 runs is fantastic.

---
 TrimSP.html        |   6 +-
 TrimSPelec.js      |  28 ++++-
 TrimSPlib.js       | 306 ++++++++++++++++++++++++++++++++-------------
 fortran/trimspNL.F |  29 ++---
 4 files changed, 259 insertions(+), 110 deletions(-)

diff --git a/TrimSP.html b/TrimSP.html
index 3bd7cc3..8229998 100644
--- a/TrimSP.html
+++ b/TrimSP.html
@@ -97,7 +97,7 @@
 	  </div>
 	  <table id="ScansTable" style="width: 100%;visibility: hidden;">
 	    <tr>
-	      <td colspan="3">Scan parameter
+	      <td colspan="1">Scan:
 		<select name="comboScan" id="comboScan" onchange="">
 		  <option value="EScan">Energy [eV]</option>
 		  <option value="SigEScan">Energy sigma [eV]</option>
@@ -107,6 +107,10 @@
 		  <option value="dScan">Thickness of layer # 1</option>
 		</select>
 	      </td>
+	      <td colspan="1">
+		Threads:
+		<input name="scanThreads" id="scanThreads" type="number" min="1" step="1" value="">
+	      </td>
 	    </tr>
 	    <tr>
 	      <td>
diff --git a/TrimSPelec.js b/TrimSPelec.js
index c524aa5..506403e 100644
--- a/TrimSPelec.js
+++ b/TrimSPelec.js
@@ -4,7 +4,8 @@
 const Plotly = require('plotly.js-dist');
 const ipcRenderer = require('electron').ipcRenderer;
 const fs = require('fs');
-const exec = require('child_process').execSync;
+const { execSync, execFile } = require('child_process');
+const os = require('os');
 const path = require('path');
 
 function getFiles(dir, filelist){
@@ -20,8 +21,27 @@ function getFiles(dir, filelist){
     return fileList;
 }
 
-function execute(command) {
-    exec(command, {stdio: 'inherit'});
+// Run one TrimSP binary asynchronously in a given working directory so the
+// scan loop can keep several local jobs in flight at once.
+function executeFileAsync(binary, args, workdir) {
+    return new Promise((resolve, reject) => {
+	execFile(binary, args, {cwd: workdir}, (err) => {
+	    if (err) {
+		reject(err);
+		return;
+	    }
+	    resolve();
+	});
+    });
+}
+
+// Use the available CPU parallelism as the default worker-pool size for local
+// scan jobs. The caller still caps this against the number of scan points.
+function getParallelJobCount() {
+    if (typeof os.availableParallelism === 'function') {
+	return Math.max(1, os.availableParallelism());
+    }
+    return Math.max(1, os.cpus().length);
 }
 
 function checkDir(directory) {
@@ -145,7 +165,7 @@ function resolveConfigWorkPath(filename, content) {
 	return null;
     }
     if (configuredWorkPath.startsWith('./')) {
-	return path.resolve(path.dirname(filename), configuredWorkPath);
+	return path.resolve(configuredWorkPath);
     }
     if (path.isAbsolute(configuredWorkPath)) {
 	return configuredWorkPath;
diff --git a/TrimSPlib.js b/TrimSPlib.js
index 0a6e2d8..2873190 100644
--- a/TrimSPlib.js
+++ b/TrimSPlib.js
@@ -1,5 +1,9 @@
 // This file contains function that are generic Javascript for TrimSP
 
+const electronExecSync = (typeof module !== 'undefined' && module.exports)
+    ? require('child_process').execSync
+    : null;
+
 // This is an object that contains all GUI ids and values
 var All = new Object();
 
@@ -3351,6 +3355,19 @@ function adjust_scans()
     var scanSeq = document.getElementById("scanSeq").checked;
     var sender = document.getElementById("scanSeq");
     var scanType = document.getElementById("scanType").value;
+    var scanThreads = document.getElementById("scanThreads");
+
+    // Default to half the available local worker slots, rounded down, while
+    // keeping at least one process. Web mode keeps the default at 1 because
+    // the CGI backend still executes scan points sequentially.
+    if (scanThreads && (scanThreads.value === "" || Number(scanThreads.value) < 1)) {
+	let availableThreads = 1;
+	if (!amIWeb() && typeof getParallelJobCount === 'function') {
+	    availableThreads = getParallelJobCount();
+	}
+	scanThreads.value = Math.max(1, Math.floor(availableThreads / 2));
+    }
+
     // consider only if scans checkbox is checked
     if (scanSeq) {
 	// Which type of scan
@@ -3453,7 +3470,7 @@ function prep_cfg(toggle) {
     valtmp = 1*(document.getElementById(strtmp).checked);
     All[strtmp]=valtmp;
     ScanSec = ScanSec + strtmp +"="+ valtmp+"\n";
-    let parScan = ["comboScan","scanType","scanFrom","scanStep","scanTo","scanList","scanListdz"];
+    let parScan = ["comboScan","scanType","scanThreads","scanFrom","scanStep","scanTo","scanList","scanListdz"];
     for (key of parScan) {
 	All[key]= document.getElementById(key).value;
 	ScanSec = ScanSec + key + "=" + All[key] + "\n";
@@ -3629,6 +3646,77 @@ function yieldToRenderer() {
     return new Promise((resolve) => setTimeout(resolve, 0));
 }
 
+// Format the live profile legend from the scan parameter currently being run.
+function formatScanLegend(scanName, scanValue) {
+    if (scanName === "E") {
+	return scanName + '=' + Number(scanValue) / 1000 + 'keV';
+    }
+    return scanName + '=' + scanValue;
+}
+
+// Merge finished per-run *.seq files into one in-memory table ordered by the
+// original scan index, not by worker completion order.
+function buildCombinedSequenceLines(results) {
+    let combinedSeqLines = [];
+    const completedResults = results.filter(Boolean).sort((a, b) => a.index - b.index);
+
+    for (const result of completedResults) {
+	if (combinedSeqLines.length === 0) {
+	    combinedSeqLines = result.seqLines.slice();
+	} else {
+	    combinedSeqLines = combinedSeqLines.concat(result.seqLines.slice(1));
+	}
+    }
+
+    return combinedSeqLines;
+}
+
+// Rebuild the live profile and fraction plots from all completed scan jobs.
+// This keeps the display responsive while preserving deterministic ordering.
+function renderCompletedScanResults(results, allValues) {
+    const completedResults = results.filter(Boolean).sort((a, b) => a.index - b.index);
+    if (completedResults.length === 0) {
+	return [];
+    }
+
+    Plotly.purge("plotRge");
+    for (const result of completedResults) {
+	Plot_xy("plotRge", result.depth, result.nmuons,
+		['Depth (nm)', 'Stopped muons (%/nm)', result.legend]);
+    }
+
+    const combinedSeqLines = buildCombinedSequenceLines(results);
+    const [combinedCols, combinedData] = parseDatContent(combinedSeqLines.join('\n'));
+    const seqInfo = getSequenceColumnInfo(combinedCols, allValues["numLayer"]);
+    const energy = combinedData[seqInfo.indexByName['Energy']];
+    const ntot = combinedData[seqInfo.indexByName['ntot']];
+    const plotSeries = [
+	{
+	    index: seqInfo.indexByName['backsc'],
+	    label: combinedCols[seqInfo.indexByName['backsc']],
+	},
+	...seqInfo.layerIndices.map((index, offset) => {
+	    let LComp = "L" + (offset + 1) + "Comp";
+	    return {
+		index: index,
+		label: allValues[LComp] || seqInfo.layerHeaders[offset],
+	    };
+	}),
+    ];
+
+    Plotly.purge("plotFrac");
+    for (const series of plotSeries) {
+	let idata = combinedData[series.index];
+	for (let i = 0; i < idata.length; i++) {
+	    idata[i] = 100 * idata[i] / ntot[i];
+	}
+	Plot_xy("plotFrac", energy, combinedData[series.index],
+		['E (keV)', 'Fraction of muons (%)', series.label]);
+    }
+
+    return combinedSeqLines;
+}
+
 async function startSequence() {
     let cmd = '';
     let trimBin = All['trimPath']+"/trimspNL";
@@ -3693,43 +3781,46 @@ async function startSequence() {
 	SValues.push(parseInt(All["valEnergy"]));
     }
 
-    // Now start the actual simulation
-    //setInterval(() => {
-    let iScan=0;
-    for (var SValue of SValues) {
-	// Update value in GUI
-	document.getElementById(ScanAttrib).value = SValue;
-	All[ScanAttrib]=SValue;
-	if ( All["SdzFlag"] == 1) {
-	    if (All["comboScan"]=="Energy") {
-		document.getElementById(ScanAttrib).value = SdzValues[iScan];
-	    } else if (All["comboScan"]=="") {
-		document.getElementById(ScanAttrib).value = SdzValues[iScan];
-	    }
-	}
-	// Update GUI progress bar
-	Progress=Math.round(Progress+90/SValues.length);
-	document.getElementById("pBar").style.width = Progress + "%";
-	document.getElementById("pBar").innerHTML = Progress + "%";
-	updateStatusBar({ run: `Running ${iScan + 1}/${SValues.length}: ${ScanName}=${SValue}` });
-	await yieldToRenderer();
-	
-	// Single run: Start 
-	let inputFileContent=CreateInpFile();
-	if (inputFileContent=="ERROR") {return(0);}
-	let FILENAME=All["fileNamePrefix"]+"_"+ScanName+SValue;
-	writeAsciiFile(All["workPath"]+"/"+FILENAME+".inp",inputFileContent);
-	if (!webOrApp) {
-	    // Prepare command and execute
-	    cmd = "cd " + All["workPath"];
-	    cmd += ";" + trimBin + ' ' + FILENAME;
-	    execute(cmd);
-	}
-	// Single run: End
+    // Freeze the initial GUI state and expand the scan into independent jobs.
+    // Each job gets its own input file name and a deterministic seed offset so
+    // the results are reproducible even when runs finish out of order.
+    const scanValues = SValues.slice();
+    const initialAll = Object.assign({}, All);
+    const initialSeed = parseInt(initialAll["ranSeed"], 10);
+    let jobs = [];
 
-	// Read stopping profiels
-	let [cols,data]= readDatFile(All["workPath"]+"/"+FILENAME+".rge");
-	// convert depth to nm and normalize stopping profile
+    for (let iScan = 0; iScan < scanValues.length; iScan++) {
+	let SValue = scanValues[iScan];
+	let jobAll = Object.assign({}, initialAll);
+	jobAll[ScanAttrib] = SValue;
+	if (!Number.isNaN(initialSeed)) {
+	    jobAll["ranSeed"] = String(initialSeed + iScan);
+	}
+	All = jobAll;
+	let inputFileContent = CreateInpFile();
+	if (inputFileContent=="ERROR") {
+	    All = initialAll;
+	    return(0);
+	}
+	let FILENAME = jobAll["fileNamePrefix"] + "_" + ScanName + SValue;
+	jobs.push({
+	    index: iScan,
+	    scanValue: SValue,
+	    filename: FILENAME,
+	    inputFileContent: inputFileContent,
+	    legend: formatScanLegend(ScanName, SValue),
+	});
+    }
+    All = initialAll;
+
+    // Results are stored by original scan index. The plots and final merged
+    // sequence file are always rebuilt from this ordered array.
+    let results = new Array(jobs.length);
+
+    // Read back the outputs of one completed run and convert them into the
+    // normalized data structures used by the live plots.
+    async function collectRunResult(job) {
+	let [cols,data]= readDatFile(All["workPath"]+"/"+job.filename+".rge");
 	let depth=data[0];
 	let nmuons=data[1];
 	let dz = (depth[1]-depth[0])/10;
@@ -3738,58 +3829,94 @@ async function startSequence() {
 	    depth[i]=depth[i]/10;
 	    nmuons[i]=nmuons[i]/norm;
 	}
-	Plot_xy("plotRge",depth,nmuons,['Depth (nm)','Stopped muons (%/nm)',ScanName+'='+Number(SValue)/1000+'keV']);		
-	// Read the sequence data
-	let [cfort,cdata]= readDatFile(All["workPath"]+"/fort.33");
-	let seqInfo = getSequenceColumnInfo(cfort, All["numLayer"]);
-	let energy = cdata[seqInfo.indexByName['Energy']];
-	let ntot = cdata[seqInfo.indexByName['ntot']];
-	let plotSeries = [
-	    {
-		index: seqInfo.indexByName['backsc'],
-		label: cfort[seqInfo.indexByName['backsc']],
-	    },
-	    ...seqInfo.layerIndices.map((index, offset) => {
-		let LComp = "L" + (offset + 1) + "Comp";
-		return {
-		    index: index,
-		    label: All[LComp] || seqInfo.layerHeaders[offset],
-		};
-	    }),
-	];
-	// Now loop over all
-	Plotly.purge("plotFrac")
-	for (let series of plotSeries) {
-	    // Normalize
-	    let idata = cdata[series.index];
-	    for (let i=0; i<idata.length; i++) {
-		idata[i]=100*idata[i]/ntot[i];
-	    }
-	    Plot_xy("plotFrac",energy,cdata[series.index],['E (keV)','Fraction of muons (%)',series.label]);
-	}
+	let seqRunFile = All["workPath"]+"/"+job.filename+".seq";
+	let seqRunLines = readAsciiFile(seqRunFile).trim().split('\n');
+	return {
+	    index: job.index,
+	    scanValue: job.scanValue,
+	    filename: job.filename,
+	    legend: job.legend,
+	    depth: depth,
+	    nmuons: nmuons,
+	    seqLines: seqRunLines,
+	};
+    }
+
+    // Refresh the progress bar and rebuild the plots from every completed job.
+    async function updateCompletedResults() {
+	const completedCount = results.filter(Boolean).length;
+	Progress = Math.round(90 * completedCount / jobs.length);
+	document.getElementById("pBar").style.width = Progress + "%";
+	document.getElementById("pBar").innerHTML = Progress + "%";
+	renderCompletedScanResults(results, initialAll);
 	await yieldToRenderer();
-	iScan++;
     }
-    
-    let data = readAsciiFile(All["workPath"]+"/"+"fort.33");
-    let LComp = "";
-    let chem_formula = "";
-    let place_holder = "";
-    let re = new RegExp(place_holder,"g");
-    for (let i=1;i<=All["numLayer"];i++) {
-	LComp = "L"+i+"Comp";
-	chem_formula = All[LComp];
-	place_holder = "impL"+i;
-	re = new RegExp(place_holder,"g");
-	data = data.replace(re, chem_formula);
+
+    // Web mode stays sequential because the CGI path still submits one job at
+    // a time. Local/Electron mode can use a bounded worker pool.
+    if (webOrApp || jobs.length <= 1 || typeof executeFileAsync !== 'function') {
+	for (const job of jobs) {
+	    writeAsciiFile(All["workPath"]+"/"+job.filename+".inp", job.inputFileContent);
+	    if (!webOrApp) {
+		cmd = "cd " + All["workPath"];
+		cmd += ";" + trimBin + ' ' + job.filename;
+		electronExecSync(cmd, {stdio: 'inherit'});
+	    }
+	    results[job.index] = await collectRunResult(job);
+	    updateStatusBar({ run: `Completed ${results.filter(Boolean).length}/${jobs.length}: ${job.legend}` });
+	    await updateCompletedResults();
+	}
+    } else {
+	// Keep a small pool of local TrimSP processes running concurrently. Jobs
+	// are inserted into the final result set by scan index as they complete.
+	let requestedThreads = parseInt(initialAll["scanThreads"], 10);
+	if (Number.isNaN(requestedThreads) || requestedThreads < 1) {
+	    requestedThreads = Math.max(1, Math.floor(getParallelJobCount() / 2));
+	}
+	const maxConcurrent = Math.min(jobs.length, requestedThreads, Math.max(1, getParallelJobCount()));
+	let nextJob = 0;
+	let activeJobs = 0;
+	let completedJobs = 0;
+
+	await new Promise((resolve, reject) => {
+	    function launchNextJobs() {
+		while (activeJobs < maxConcurrent && nextJob < jobs.length) {
+		    const job = jobs[nextJob];
+		    nextJob += 1;
+		    activeJobs += 1;
+		    writeAsciiFile(All["workPath"]+"/"+job.filename+".inp", job.inputFileContent);
+		    updateStatusBar({
+			run: `Running ${completedJobs}/${jobs.length} completed, ${activeJobs} active`
+		    });
+
+		    executeFileAsync(trimBin, [job.filename], All["workPath"])
+			.then(() => collectRunResult(job))
+			.then(async (result) => {
+			    results[result.index] = result;
+			    activeJobs -= 1;
+			    completedJobs += 1;
+			    updateStatusBar({
+				run: `Completed ${completedJobs}/${jobs.length}, ${activeJobs} active`
+			    });
+			    await updateCompletedResults();
+			    if (completedJobs === jobs.length) {
+				resolve();
+				return;
+			    }
+			    launchNextJobs();
+			})
+			.catch(reject);
+		}
+	    }
+
+	    launchNextJobs();
+	});
     }
+
+    let combinedSeqLines = buildCombinedSequenceLines(results);
     let seq_file = All["workPath"]+"/"+All["fileNamePrefix"]+"_Seq_Results.dat";
-    writeAsciiFile(seq_file, data);
-    
-    // Remove redundant files and change the name fort.33
-    if (!webOrApp) {
-	cmd="cd " + All["workPath"] + "; mv -f fort.33 " + seq_file;
-	execute(cmd);
+    if (combinedSeqLines.length > 0) {
+	writeAsciiFile(seq_file, combinedSeqLines.join('\n') + '\n');
     }
     Progress=100;
     document.getElementById("pBar").style.width = Progress + "%";
@@ -3803,9 +3930,12 @@ function readDatFile(filename) {
     // Read column data file and return
     // cols - labels of columns
     // data - 2D array with the columns
-    let lines = "";
-    lines = readAsciiFile(filename);
-    
+    let lines = readAsciiFile(filename);
+    return parseDatContent(lines);
+}
+
+function parseDatContent(lines) {
+    // Parse whitespace-separated column data from an in-memory string.
     let data = [];
     let i=0; // line counter
     for (let line of lines.trim().split('\n')) {
diff --git a/fortran/trimspNL.F b/fortran/trimspNL.F
index 39b15e7..164241d 100644
--- a/fortran/trimspNL.F
+++ b/fortran/trimspNL.F
@@ -71,7 +71,7 @@ C     Maximum number of elements in each layer, was limited to 5.
       PARAMETER (MAXNLm15=(MAXNL-1)*MAXEL)
       PARAMETER (TRIMSP_VERSION='1.3.2')
       LOGICAL   TEST(64),TESTR(2000),TEST1(2000)
-      LOGICAL   EQUAL,FORT33
+      LOGICAL   EQUAL
       INTEGER*4 ISRCHFGT,ISRCHFGE,ILLZ
       INTEGER   N,L,LL,NH,NUM,KK
       INTEGER   I,J,IV
@@ -284,8 +284,8 @@ c     for part. reflec. coeff. calculation
 C     CHARACTER Variables
       CHARACTER*18 DPOT,DPOTR,DKDEE1,DKDEE2
       CHARACTER*60 filein,fileout
-      CHARACTER*64 innam,outnam,rgenam,errnam
-      CHARACTER*4 inext,outext,rgeext,errext
+      CHARACTER*64 innam,outnam,rgenam,errnam,seqnam
+      CHARACTER*4 inext,outext,rgeext,errext,seqext
       CHARACTER errcom*72
       CHARACTER month_start*4,month_stop*4,day_start*2,day_stop*2
       CHARACTER year_start*4,year_stop*4,hour_start*2,hour_stop*2
@@ -298,7 +298,7 @@ C     CHARACTER Variables
       DATA  PI/3.14159265358979D0/, ICW/100/, E2/14.399651D0/
       DATA  AB/0.52917725D0/, FP/0.885341377D0/, AN/0.60221367D0/
       DATA  inext/'.inp'/,outext/'.out'/,rgeext/'.rge'/
-      DATA  errext/'.err'/
+      DATA  errext/'.err'/,seqext/'.seq'/
 
       DATA  ET/0.D0/,PLST/0.D0/,PL2ST/0.D0/,PL3ST/0.D0/
       DATA  PL4ST/0.D0/,PL5ST/0.D0/,PL6ST/0.D0/
@@ -431,6 +431,7 @@ C     Require an explicit run basename for input/output naming.
       outnam=trim(fileout)//outext
       rgenam=trim(fileout)//rgeext
       errnam=trim(fileout)//errext
+      seqnam=trim(fileout)//seqext
 
       write (*,*) innam
 
@@ -3465,7 +3466,7 @@ C
  1800 CONTINUE
 
 C     
-C     The file fort.33 is created here
+C     The per-run summary file (*.seq) is created here
 C     This is the summary file writing procedure
 C
  7802 FORMAT(6x,'Energy',4x,'SigmaE',5x,'Alpha',2x,'SigAlpha',4x,'ntot',
@@ -3477,23 +3478,17 @@ C     &     I0,3x))
  7801 FORMAT(F12.2,3(1x,F9.2),1x,6(I7,1x),6(E12.4),2(E12.4),999(I7
      &     ,1x))
 
-      inquire(FILE='fort.33',EXIST=FORT33)
-C     fort.33 is the compact per-run summary consumed by the GUI for scan
-C     plots (fractions stopped in each layer, backscattering, transmission,
-C     mean depth, etc.). The header line is written once, then each run
-C     appends a single summary row.
-      if (.not.FORT33) then
-         open(33)
-         WRITE(33,7802)  (chem(k),k=1,NLayers)
-      else
-         open(33,access='append')
-      endif
+C     The *.seq file is the compact per-run summary consumed by the GUI for
+C     scan plots (fractions stopped in each layer, backscattering,
+C     transmission, mean depth, etc.). Each run writes its own summary file.
+      OPEN(UNIT=33,FILE=seqnam,STATUS='replace')
+      WRITE(33,7802)  (chem(k),k=1,NLayers)
       WRITE(33,7801)E0keV,EsigkeV,ALPHA,ALPHASIG,NH,IIM,IB,IT,tryE
      &     ,negE,FIX0,SIGMAX,FIB0,SIGMAB,FIT0,SIGMAT,epsilon,prcoeff
      &     ,(number_in_layer(k),k=1,NLayers)
       CLOSE(33)
 C
-C     End of file fort.33
+C     End of per-run summary file
       
 C
 C     TOP AND FRONT LINES FOR MATRICES