Synthetic turbulence: draw modes from an exact integer PRNG so seeds are independent

[sbryngelson]

Description

Different synth_seed values produce mostly the same synthetic-turbulence forcing, so seeded ensembles are not independent.

Root cause. s_initialize_body_forces_module draws the wave vectors and phases with s_prng/modmul (src/common/m_helper.fpp). modmul rounds the LCG state to 5 decimals every step (decimal_trim = 1e5). That makes the map many-to-one: every seed falls into the same short cycle, and streams from different seeds merge. Here is a Python reproduction of s_prng (same constants and rounding):

seed 1:     enters cycle after 385 draws, cycle length 149, cycle min 161061
seed 2:     enters cycle after  50 draws, cycle length 149, cycle min 161061
seed 3:     enters cycle after 448 draws, cycle length 149, cycle min 161061
seed 12345: enters cycle after 132 draws, cycle length 149, cycle min 161061

The same holds for seeds 4, 5, 10 and 100. In 3D each mode takes 5 draws, so 48 modes use 240 draws:

s_prng    : first draw seeds 1,2,3 = [0.946, 0.947, 0.949]; max shared values among first 240 draws of any 2 seeds = 128
splitmix32: first draw seeds 1,2,3 = [0.588, 0.704, 0.928]; max shared values among first 240 draws of any 2 seeds = 0
splitmix32 seed 1, 1e6 draws: mean 0.5003, var 0.0832 (1/12 = 0.0833), distinct 1000000

Seeds tested: 1, 2, 3, 4, 5, 10, 100, 12345.

Observed impact. We ran 3D forced flow past a flapping wing on Frontier with seeds 1, 2 and 3 as an "independent" ensemble. The upstream probe's w(t) correlated at r = 0.44-0.76 between seeds over t = 3-5, and at r = 0.73-0.85 over t = 4-12. Lift during the first wingbeat had a seed-to-seed standard deviation of 0.003, against an excursion of 1.7.

Fix.

• Add s_prng_splitmix32 to m_helper: a 32-bit SplitMix-style generator, i.e. a Weyl sequence with step 0x9E3779B9 hashed by the MurmurHash3 fmix32 finalizer. It has period 2^32 and is a bijective hash, so there is no cycle collapse.
• All arithmetic is exact 64-bit integer (selected_int_kind(18), as elsewhere in m_helper). The 32-bit multiply is split over 16-bit halves (f_mulmod32), so no intermediate exceeds 2^49. There is no signed overflow and no floating-point rounding, so the stream is identical on every compiler. I checked this with a standalone program: gfortran and CCE give bit-identical output to each other and to a Python model.
• Use it for the synthetic-turbulence mode draws. seed becomes a 64-bit integer initialized from synth_seed, which the generator reduces mod 2^32. The generator runs only on rank 0 at init, on the host, so there are no GPU implications.

s_prng itself is unchanged. Its only other caller is s_generate_random_perturbation in src/pre_process/m_perturbation.fpp (the mixing-layer perturbation, test mixlayer_perturb). That caller reseeds per mode and y-location and takes only 5 draws per seed. Changing s_prng would silently change those initial conditions and goldens, so this PR leaves it alone. The rounding can still map different (k, y) inputs onto correlated draws there. Maintainers may want to switch that caller too, as a separate change.

Impact on results / goldens

• The realized forcing changes for every existing synth_seed, so runs with synthetic turbulence will not reproduce their pre-PR forcing. Nothing else changes.
• Golden CB853530 (3D -> synthetic_turbulence, synth_seed = 1) is regenerated in this PR (commit 80b58e1). It was generated on a Frontier CPU compute node with GNU Fortran 12.3.0 (gcc-native/12.3) in Debug mode without MPI, via ./mfc.sh test --generate --only CB853530 --no-gpu --no-mpi --debug. This matches the configuration recorded in the old golden's metadata (GNU 12.2.0, Debug, no MPI) except for the GCC minor version.
  • Only tests/CB853530/ changed. The new golden has the same 12 entries and sizes as the old one. The t = 0 fields are identical; the step-3 fields differ by at most about 1e-3, the size of the forcing increment itself, and all values are finite.
  • Against the new golden, ./mfc.sh test --only CB853530 passes.
  • D13FDB23 (3D -> mixlayer_perturb, the other s_prng user) still passes against its unchanged golden.
• examples/2D_synthetic_turbulence is already in the skip list in cases.py, so its golden is not checked.

Verification

• The cycle and coalescence numbers above come from a Python model of s_prng and of the new generator.
• s_prng_splitmix32 in standalone Fortran matches the Python model bit-for-bit (first 3 draws for seeds 1-3, and draw 10^6 for seed 12345) under both gfortran and CCE 19.
• simulation builds on CPU with CCE 19 and with GNU 12.3. Apart from the 3-step regression tests above, I did not run this generator in a simulation. In our production runs we confirmed independence with a local, opt-in workaround that used the intrinsic random_number; I am not proposing that here because it is not reproducible across compilers.

Contribution Policy

We do not accept pull requests generated primarily by AI without genuine understanding or real-world usage context.

All contributions are expected to demonstrate:

• A clear understanding of the codebase
• Alignment with product direction
• Thoughtful reasoning behind changes
• Evidence of real-world usage or hands-on experience with the problem

If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.

---

Acknowledgement

• I confirm this PR meets the above expectations and reflects my own understanding and real-world context.

---

This PR was prepared with the assistance of an AI tool (Claude Code). We found the problem while analyzing a seeded ensemble of production runs on Frontier.

PR template credit: junegunn

[Copilot]

Warning

Copilot couldn't run its full agentic review because it didn't start before the timeout. Make sure your repository has a runner available, or add a copilot-code-review.yml file specifying one with the runs-on attribute. See the docs for more details.

Copilot review overview

Review effort: Lite
Findings: 4 · 2 · 2

Open (8)

• This change alters the realized synthetic-turbulence forcing for existing synth_seed values, so… · New
• This change alters the realized synthetic-turbulence forcing for existing synth_seed values, so… · New
• This change alters the realized synthetic-turbulence forcing for existing synth_seed values, so… · New
• This change alters the realized synthetic-turbulence forcing for existing synth_seed values, so… · New
• Using integer(kind=8) relies on compiler-specific kind numbering. For portability (and to match… · New
• Using integer(kind=8) relies on compiler-specific kind numbering. For portability (and to match… · New
• The routine returns var = real(z, wp)/4294967296._wp, which yields values in [0, 1) (1.0 is not… · New
• This comment line is very long and may exceed typical free-form Fortran line-length limits used in… · New

What changed in this PR

This PR fixes synthetic-turbulence seeding by replacing the current rounded LCG-based draws with an exact 32-bit integer SplitMix-style PRNG to prevent seed stream collapse and improve ensemble independence while keeping compiler-independent reproducibility.

Changes:

• Switch synthetic-turbulence mode initialization to use s_prng_splitmix32 and a 64-bit integer PRNG state masked to 32 bits.
• Add s_prng_splitmix32 (and helper f_mulmod32) to m_helper implementing exact integer arithmetic with a MurmurHash3-style finalizer.

File	Description
src/​simulation/​m_body_forces.fpp	Uses the new SplitMix-style PRNG for synthetic turbulence mode vector/phase draws and updates seed handling.
src/​common/​m_helper.fpp	Introduces s_prng_splitmix32 (SplitMix/Weyl + fmix32) and f_mulmod32 for exact 32-bit modular multiplication.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[sbryngelson]

Done in 4d6c129. I used selected_int_kind(18), the convention m_helper already uses for 64-bit integers (int64_kind in double_factorial/factorial); src/ does not use iso_fortran_env for this. The literals are now _int64_kind, and the caller in m_body_forces declares seed with the same kind. The generator also reduces the state mod 2^32 itself, so the caller no longer needs a kind-specific mask literal. The output stream is bit-identical to before: I checked against a Python model with gfortran and CCE, including negative seeds.

[sbryngelson]

Done in 4d6c129. I used selected_int_kind(18), the convention m_helper already uses for 64-bit integers (int64_kind in double_factorial/factorial); src/ does not use iso_fortran_env for this. The literals are now _int64_kind, and the caller in m_body_forces declares seed with the same kind. The generator also reduces the state mod 2^32 itself, so the caller no longer needs a kind-specific mask literal. The output stream is bit-identical to before: I checked against a Python model with gfortran and CCE, including negative seeds.

[sbryngelson]

Fixed in 4d6c129: the docstring now says [0, 1). One caveat: with wp in single precision, real(2^32-1, sp)/2^32 rounds to exactly 1.0, so the comment says 1 is reachable only through that rounding.

[sbryngelson]

I wrapped it in 4d6c129, but ./mfc.sh format (ffmt) reflows comment blocks to fill the repo's 132-column limit, so it rejoins manual wraps. The resulting lines are at most 127 columns, within the limit and like the surrounding code, and the precheck formatting check passes.

[github-actions]

Lines of Code

File	Lines	Diff
src/common/m_helper.fpp	507	+17

Directory	Lines	Diff
common	10443	+17
total	47025	+17

[codecov]

Codecov Report

❌ Patch coverage is 94.11765% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 62.66%. Comparing base (3dc5b2f) to head (4d6c129).

Files with missing lines	Patch %	Lines
src/simulation/m_body_forces.fpp	85.71%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1946      +/-   ##
==========================================
+ Coverage   62.64%   62.66%   +0.01%     
==========================================
  Files          86       86              
  Lines       22425    22435      +10     
  Branches     3325     3325              
==========================================
+ Hits        14048    14058      +10     
  Misses       6119     6119              
  Partials     2258     2258              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.