cogen: decorrelate per-ray KISS seeds across seeds, MPI ranks and batches (fixes #2807)

.github/workflows/mcstas-basictest.yml

@@ -404,7 +404,7 @@ jobs:
            mcstas${EXESUFFIX} --version
            mctest${SUFFIX} --verbose ${NOPLOTS} --instr=Test_RNG_rand01 --testdir=run_rngtest --suffix=${{ matrix.cogen }}_${{ matrix.method }} --uid=${MYWORKERID}
            export SUM=`find run_rngtest -name rngout.dat | xargs -n1 grep -v \# | ${MD5SUM} | cut -f1 -d\ `
-           export EXPECTED="f192ce4609e2225bf9d42ce9c5fa5a86"
+           export EXPECTED="6c203d69205b0656c1391555e6068c93"
            if [ "${EXPECTED}" == "${SUM}" ];
            then
               echo RNG test success!

.github/workflows/mcxtrace-basictest.yml

@@ -400,7 +400,7 @@ jobs:
            mcxtrace${EXESUFFIX} --version
            mxtest${SUFFIX} --verbose ${NOPLOTS} --instr=Test_RNG_rand01 --testdir=run_rngtest --suffix=${{ matrix.cogen }}_${{ matrix.method }} --uid=${MYWORKERID}
            export SUM=`find run_rngtest -name rngout.dat | xargs -n1 grep -v \# | ${MD5SUM} | cut -f1 -d\ `
-           export EXPECTED="f192ce4609e2225bf9d42ce9c5fa5a86"
+           export EXPECTED="6c203d69205b0656c1391555e6068c93"
            if [ "${EXPECTED}" == "${SUM}" ];
            then
               echo RNG test success!

common/lib/share/mccode-r.c

@@ -4254,7 +4254,6 @@ long sort_absorb_last(_class_particle* particles, _class_particle* pbuffer, long
   // copy non-absorbed block
   #pragma acc parallel loop present(particles[0:buffer_len])
   for (long tidx = 0; tidx < accumlen; tidx++) { // tidx: thread index
-    randstate_t randstate[7];
     _class_particle sourcebuffer;
     _class_particle targetbuffer;
     // assign reduced weight to all particles
@@ -4266,13 +4265,8 @@ long sort_absorb_last(_class_particle* particles, _class_particle* pbuffer, long
       // buffer full particle struct
       targetbuffer = particles[tidx];
       // reassign previous randstate
-      targetbuffer.randstate[0] = sourcebuffer.randstate[0];
-      targetbuffer.randstate[1] = sourcebuffer.randstate[1];
-      targetbuffer.randstate[2] = sourcebuffer.randstate[2];
-      targetbuffer.randstate[3] = sourcebuffer.randstate[3];
-      targetbuffer.randstate[4] = sourcebuffer.randstate[4];
-      targetbuffer.randstate[5] = sourcebuffer.randstate[5];
-      targetbuffer.randstate[6] = sourcebuffer.randstate[6];
+      for (int ridx = 0; ridx < RANDSTATE_LEN; ridx++)
+        targetbuffer.randstate[ridx] = sourcebuffer.randstate[ridx];
       // apply
       particles[bidx*accumlen + tidx] = targetbuffer;
     }

docs/GRAMMAR/ADR-records/ADR_20261018_PER_RAY_SEEDING.md

@@ -0,0 +1,74 @@
+# Per-ray RNG seeding and `_uid` in the generated ray-trace loops
+
+## Status
+
+*Proposed* and *implemented* in mccode-dev/McCode#2808 (fixes
+mccode-dev/McCode#2807). Needs a matching change in mccode-antlr (mccode-dev/mccode-antlr#392), whose
+`translators/c_raytrace.py` holds its own copy of these lines.
+
+## Context
+
+`raytrace_all()` and `raytrace_all_funnel()` (written by `cogen.c.in`)
+seed each ray with `srandom(_hash((pidx+1)*(seed+1)))`, and then do
+`seed += gpu_innerloop` after every batch. Because the key is a product:
+
+* runs at different seeds share rays: ray `p` at seed 3 is ray `2p+1`
+  at seed 1 (50% of rays shared);
+* the ranks of one MPI run (`mcseed + rank`) share rays: with
+  `-np 8 --seed 1` only ~55% of rays are distinct;
+* with 32-bit `unsigned long` (MSVC) the key wraps at 2^32; seed
+  2147483647 gives just two different rays.
+
+Also, `_uid` was `pidx` (+ `rank*ncount` with MPI), so it repeated in
+every GPU batch. The rank stride `ncount` could also be smaller than
+the `loops*gpu_innerloop` rays a rank actually traces.
+
+## Decision
+
+In both loops (KISS RNG, `RNG_ALG=2`, the default):
+
+```c
+_particle->_uid = cloop*gpu_innerloop + pidx;
+_particle->_uid += mpi_node_rank * loops*gpu_innerloop;          /* USE_MPI */
+srandom(_hash(_hash(seed) + cloop*gpu_innerloop + pidx));
+```
+
+and drop the per-batch `seed += gpu_innerloop`. The global ray index
+replaces it. Keeping it would make seed S batch 1 equal to seed
+S+gpu_innerloop batch 0.
+
+`_hash` is the existing 64-bit bijection (`randstate_t` is `uint64_t`
+for KISS) and is already `#pragma acc routine seq`. No runtime change and
+no new helper: a SplitMix64 mixer was tried and gave the same results.
+
+## Consequences
+
+* **Numerical results change** for every KISS run at a given `--seed`.
+  Outputs are not bit-identical to earlier releases. The CI
+  `Test_RNG_rand01` checksum was updated to
+  `6c203d69205b0656c1391555e6068c93`. Mersenne Twister builds
+  (`RNG_ALG=1`) ignore the per-ray seed and are unaffected.
+* Different seeds and different MPI ranks no longer share rays. MSVC
+  (32-bit `long`) behaves like the other platforms.
+* On CPU (`loops=1`, `gpu_innerloop=ncount`) `_uid` values are unchanged.
+  On GPU they are unique across batches.
+* The particle's `randstate` is now sized by `RANDSTATE_LEN` (was a hard-coded
+  `[7]`), and `sort_absorb_last()` (FUNNEL + SPLIT) copies `RANDSTATE_LEN` words
+  instead of 7 hand-written ones. No change for KISS (7); the MT struct shrinks to
+  1 word. This lets an RNG with another state length be added without
+  out-of-bounds copies.
+* mccode-antlr must port the same lines by hand (mccode-dev/mccode-antlr#392). Until it does, the two
+  code generators give different results for the same seed.
+
+## Behaviour
+
+Measured on macOS arm64, clang, CPU (rays identified by their first draws):
+
+| case | before | after |
+|---|---|---|
+| seeds 1 and 3, rays shared | 50% | 0 |
+| seed 2147483647, 32-bit, distinct rays | 2 | all |
+| `mpirun -np 8 --seed 1`, distinct rays | ~55% | 100% |
+| same seed twice | identical | identical |
+
+Not tested locally: MSVC, Linux, nvc/OpenACC (batched GPU path).

mccode/src/cogen.c.in

@@ -2140,16 +2140,15 @@ int cogen_raytrace(struct instr_def *instr)
   coutf("    for (unsigned long pidx=0 ; pidx < gpu_innerloop ; pidx++) {");
   coutf("      _class_particle particleN = mcgenstate(); // initial particle");
   coutf("      _class_particle* _particle = &particleN;");
-  coutf("      particleN._uid = pidx;");
+  coutf("      particleN._uid = cloop*gpu_innerloop + pidx;");
   coutf("      #ifdef USE_MPI");
-  coutf("      particleN._uid += mpi_node_rank * ncount; ");
+  coutf("      particleN._uid += mpi_node_rank * loops*gpu_innerloop;");
   coutf("      #endif");
   coutf("");
-  coutf("      srandom(_hash((pidx+1)*(seed+1)));");
+  coutf("      srandom(_hash(_hash(seed) + cloop*gpu_innerloop + pidx));");
   coutf("");
   coutf("      raytrace(_particle);");
   coutf("    } /* inner for */");
-  coutf("    seed = seed+gpu_innerloop;");
   coutf("  } /* CPU for */");
   coutf("  /* if on GPU, printf has been globally nullified, re-enable here */");
   cout("     #ifdef OPENACC");
@@ -2279,11 +2278,11 @@ int cogen_rt_funnel(struct instr_def *instr)
   coutf("      // generate particle state, set loop index and seed");
   coutf("      particles[pidx] = mcgenstate();");
   coutf("      _class_particle* _particle = particles + pidx;");
-  coutf("      _particle->_uid = pidx;");
+  coutf("      _particle->_uid = cloop*gpu_innerloop + pidx;");
   coutf("      #ifdef USE_MPI");
-  coutf("      _particle->_uid += mpi_node_rank * ncount; ");
+  coutf("      _particle->_uid += mpi_node_rank * loops*gpu_innerloop;");
   coutf("      #endif");
-  coutf("      srandom(_hash((pidx+1)*(seed+1))); // _particle->state usage built into srandom macro");
+  coutf("      srandom(_hash(_hash(seed) + cloop*gpu_innerloop + pidx)); // _particle->state usage built into srandom macro");
   coutf("    }");
   cout( "");
 
@@ -2394,8 +2393,6 @@ int cogen_rt_funnel(struct instr_def *instr)
   cout( "");
   list_iterate_end(liter);
 
-  coutf("    // jump to next viable seed");
-  coutf("    seed = seed + gpu_innerloop;");
   coutf("  } // outer loop / particle batches");
   cout( "");
   coutf("  free(particles);");
@@ -2567,8 +2564,7 @@ cogen_header(struct instr_def *instr, char *output_name)
   cout("  double _mctmp_a; /* temp a */");
   cout("  double _mctmp_b; /* temp b */");
   cout("  double _mctmp_c; /* temp c */");
-  //cout("  randstate_t randstate[RANDSTATE_LEN]");
-  cout("  randstate_t randstate[7];"); // for the KISS generator
+  cout("  randstate_t randstate[RANDSTATE_LEN];");
   cout("  double t, p;     /* time, event weight */");
   cout("  long long _uid;  /* Unique event ID */");
   cout("  long _index;     /* component index where to send this event */");

mcstas-comps/examples/FZ_Juelich/FZJ_BenchmarkSfin2/FZJ_BenchmarkSfin2.instr

@@ -15,7 +15,7 @@
 * Test instrument for the H Frielinghaus SANS_benchmark2 component. 
 *
 * The below test uses defaults of the instrument (SANS_benchmark2 modnum=11)
-* %Example: -y Detector: PSDnear_I=7400
+* %Example: -y Detector: PSDnear_I=8.5e3
 *
 * Many sample-modes are available, from the SANS_benchmark2 component:
 * The component includes 19 sample-model-settings:

mcstas-comps/examples/Mantid/templateSANS_Mantid/templateSANS_Mantid.instr

@@ -26,8 +26,8 @@
 * 2) Generate an IDF using mcdisplay templateSANS_Mantid --format=Mantid -n0
 * 3) mcrun templateSANS_Mantid -n1e6 --format=NeXus
 *
-* %Example: lambda=6 Detector: PSDrad_I=105000000
-* %Example: lambda=6 --format=NeXus Detector: PSDrad_I=105000000
+* %Example: lambda=6 Detector: PSDrad_I=9.6e7
+* %Example: lambda=6 --format=NeXus Detector: PSDrad_I=9.6e7
 *
 * %Parameters
 * lambda: [AA]          Mean wavelength of neutrons
@@ -78,7 +78,7 @@ COMPONENT LdetectorPRE = L_monitor(
 AT (0,0,0.05) RELATIVE coll2
 
 
-SPLIT COMPONENT sampleMantid = Sans_spheres(
+SPLIT 100 COMPONENT sampleMantid = Sans_spheres(
   R=r, Phi=PHI, Delta_rho=Delta_Rho, sigma_abs=sigma_abs,
   xwidth=0.01, yheight=0.01, zdepth=0.005)
    AT (0,0,0.2) RELATIVE coll2

mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr

@@ -52,11 +52,11 @@
 * %Example: h=1.2 l=1 dE=5 sample_choice=0 primitive_cell=1 Detector: detector_I=1.46e-11
 * %Example: h=1.15 l=0 dE=3 T=2 sample_choice=2 Detector: detector_I=8.91e-10
 * %Example: h=1.15 l=0 dE=3 T=2 sample_choice=3 Detector: detector_I=9.41e-10
-* %Scan: mcrun Test_Dispersion_relation_TAS.instr -M -N 5,5 -n 2e5 h=0.6,1.4 l=0.6,1.4 dE=5 sample_choice=0 Detector: detector_I={
-*   0.0,0.0,0.0,0.0,0.0,0.0,1.2167e-12,3.52579e-11,
-*   2.26406e-12,0.0,0.0,1.17383e-11,0.0,1.39525e-11,0.0,0.0,
-*   6.10596e-12,1.00579e-11,8.66923e-12,0.0,0.0,0.0,0.0,0.0,
-*   0.0 }
+* %Scan: mcrun Test_Dispersion_relation_TAS.instr -M -N 5,5 -n 2e5 h=0.8,1.2 l=1.1,1.3 dE=5 sample_choice=0 Detector: detector_I={
+*   8.964e-12,1.086e-11,8.771e-12,1.261e-11,1.151e-11,7.748e-12,1.057e-11,1.137e-11,
+*   1.199e-11,1.014e-11,6.126e-12,8.926e-12,1.006e-11,1.007e-11,8.424e-12,3.55e-12,
+*   7.184e-12,8.291e-12,8.243e-12,5.35e-12,9.55e-13,4.109e-12,5.612e-12,4.924e-12,
+*   1.74e-12 }
 *
 * %Parameters
 * h: [r.l.u.]       Scattering vector component along a* (sample x axis)

mcstas-comps/examples/Tests_samples/Test_Magnon_bcc_TAS/Test_Magnon_bcc_TAS.instr

@@ -12,7 +12,7 @@
 * %D
 * Generic TAS instrument for test of samples with dispersions. Modeled over the RITA-2 instrument at PSI (one analyzer only). The instrument is able to position in q-space at q=(h 0 0) and fixed Ei and Ef. This can be used to longitudinal constant-E scans or constant-q scans. In addition, transverse constant-E scans can be made by scanning the sample orientation A3.
 * 
-* %Example: Test_Magnon_bcc_TAS.instr hw=0.5 Detector: e_monitor1_I=460000
+* %Example: Test_Magnon_bcc_TAS.instr hw=0.5 Detector: e_monitor1_I=4.97e5
 * %Scan: mcrun Test_Magnon_bcc_TAS.instr hw=0,1 -N21 -n1e6 Detector: e_monitor_I={
 *   5541440.0,5681320.0,4956300.0,4947280.0,4886770.0,4370310.0,3596940.0,3900100.0,
 *   3656700.0,3712310.0,3128430.0,3254960.0,2708640.0,2961760.0,2536780.0,2524220.0,
@@ -149,7 +149,7 @@ COMPONENT l_monitor = L_monitor(
     restore_neutron=1)
 AT (0, 0, 1.2) RELATIVE arm
 
-COMPONENT magnon_simple = Magnon_bcc(
+SPLIT 20 COMPONENT magnon_simple = Magnon_bcc(
     radius=sample_r,
     yheight=sample_h,
     sigma_abs=0, 

mcstas-comps/examples/Tests_samples/Test_SANS/Test_SANS.instr

@@ -28,7 +28,7 @@
 * SANSNanodiscsWithTagsFast SAMPLE=10
 * SANSPDBFast               SAMPLE=11
 *
-* %Example: SAMPLE=0 Detector: PSDMonitor_I=2.529e-09
+* %Example: SAMPLE=0 Detector: PSDMonitor_I=2.9e-09
 * %Example: SAMPLE=1 Detector: PSDMonitor_I=2.486e-10
 * %Example: SAMPLE=2 Detector: PSDMonitor_I=2.297e-10
 * %Example: SAMPLE=3 Detector: PSDMonitor_I=8.720e-11

mcstas-comps/examples/Tests_samples/Test_TasReso/Test_TasReso.instr

@@ -30,7 +30,7 @@
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see http://www.gnu.org/licenses/.
 *
-* %Example: src_lam=4.5 Detector: emon_I=887.137
+* %Example: src_lam=4.5 Detector: emon_I=960
 *
 * %Parameters 
 * src_lam:     [AA] Source mean wavelength
@@ -277,7 +277,7 @@ ROTATED (0, 0, 0) RELATIVE mono_arm2
 /* ----------------------------------------------------------------------------- */
 /* Samples */
 
-COMPONENT sample_arm1 = Arm()
+SPLIT 50 COMPONENT sample_arm1 = Arm()
 AT(0, 0, dist_mono_sample) RELATIVE mono_arm2
 ROTATED (0, 0, 0) RELATIVE mono_arm2
 

mcxtrace-comps/examples/Tests_samples/Test_SAXS/Test_SAXS.instr

@@ -35,7 +35,7 @@
 * %Example: SAMPLE=4 Detector: PSDMonitor_I=1.16433e-10
 * %Example: SAMPLE=5 Ncount=1e4 Detector: PSDMonitor_I=1.2001e-14
 * %Example: SAMPLE=6 Ncount=1e4 Detector: PSDMonitor_I=1.2400e-14
-* %Example: SAMPLE=7 Ncount=1e4 Detector: PSDMonitor_I=7.07865e-14
+* %Example: SAMPLE=7 Ncount=1e4 Detector: PSDMonitor_I=5.54e-14
 * %Example: SAMPLE=9 Detector: PSDMonitor_I=1.11073e-14
 * %Example: SAMPLE=10 Detector: PSDMonitor_I=1.02884e-14
 * %Example: SAMPLE=11 Detector: PSDMonitor_I=5.47707e-15