diff --git a/.github/workflows/mcstas-basictest.yml b/.github/workflows/mcstas-basictest.yml index 195fce8092..985b609f41 100644 --- a/.github/workflows/mcstas-basictest.yml +++ b/.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! diff --git a/.github/workflows/mcxtrace-basictest.yml b/.github/workflows/mcxtrace-basictest.yml index f760bc0e68..139e9e8ef6 100644 --- a/.github/workflows/mcxtrace-basictest.yml +++ b/.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! diff --git a/common/lib/share/mccode-r.c b/common/lib/share/mccode-r.c index c00ef93ab0..7e1903c4a6 100644 --- a/common/lib/share/mccode-r.c +++ b/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; } diff --git a/docs/GRAMMAR/ADR-records/ADR_20261018_PER_RAY_SEEDING.md b/docs/GRAMMAR/ADR-records/ADR_20261018_PER_RAY_SEEDING.md new file mode 100644 index 0000000000..d2bfe5ca47 --- /dev/null +++ b/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). diff --git a/mccode/src/cogen.c.in b/mccode/src/cogen.c.in index a322c31059..ca24eff9fc 100644 --- a/mccode/src/cogen.c.in +++ b/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 */"); diff --git a/mcstas-comps/examples/FZ_Juelich/FZJ_BenchmarkSfin2/FZJ_BenchmarkSfin2.instr b/mcstas-comps/examples/FZ_Juelich/FZJ_BenchmarkSfin2/FZJ_BenchmarkSfin2.instr index f898bb3693..3eacd2cd50 100644 --- a/mcstas-comps/examples/FZ_Juelich/FZJ_BenchmarkSfin2/FZJ_BenchmarkSfin2.instr +++ b/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: diff --git a/mcstas-comps/examples/Mantid/templateSANS_Mantid/templateSANS_Mantid.instr b/mcstas-comps/examples/Mantid/templateSANS_Mantid/templateSANS_Mantid.instr index 004c87b038..b077d577ce 100644 --- a/mcstas-comps/examples/Mantid/templateSANS_Mantid/templateSANS_Mantid.instr +++ b/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 diff --git a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr index 46b4fe57ad..93d4dea149 100644 --- a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr +++ b/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) diff --git a/mcstas-comps/examples/Tests_samples/Test_Magnon_bcc_TAS/Test_Magnon_bcc_TAS.instr b/mcstas-comps/examples/Tests_samples/Test_Magnon_bcc_TAS/Test_Magnon_bcc_TAS.instr index 5e079410c0..769b0d5801 100644 --- a/mcstas-comps/examples/Tests_samples/Test_Magnon_bcc_TAS/Test_Magnon_bcc_TAS.instr +++ b/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, diff --git a/mcstas-comps/examples/Tests_samples/Test_SANS/Test_SANS.instr b/mcstas-comps/examples/Tests_samples/Test_SANS/Test_SANS.instr index 985fb7a9ae..e842e74b8f 100644 --- a/mcstas-comps/examples/Tests_samples/Test_SANS/Test_SANS.instr +++ b/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 diff --git a/mcstas-comps/examples/Tests_samples/Test_TasReso/Test_TasReso.instr b/mcstas-comps/examples/Tests_samples/Test_TasReso/Test_TasReso.instr index 0b17bee525..dfe8347999 100644 --- a/mcstas-comps/examples/Tests_samples/Test_TasReso/Test_TasReso.instr +++ b/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 diff --git a/mcxtrace-comps/examples/Tests_samples/Test_SAXS/Test_SAXS.instr b/mcxtrace-comps/examples/Tests_samples/Test_SAXS/Test_SAXS.instr index 7f912f0a74..0b8a89803b 100644 --- a/mcxtrace-comps/examples/Tests_samples/Test_SAXS/Test_SAXS.instr +++ b/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