Synthetic turbulence

[danieljvickers]

Description

This branch includes the ability to inject synthetic turbulence in a Gaussian forcing region as described in

Tangermann, E. & Klein, M. (2020). "Controlled Synthetic Freestream Turbulence Intensity Introduced by a Local Volume Force." Fluids 5(3), 130. https://doi.org/10.3390/fluids5030130

Type of change (delete unused ones)

• New feature

Testing

I ran 2D and 3D synthetically injected turbulence over an airfoil IB.

I still need to do an analysis of the energy distribution to validate the model and show that it was implemented correctly. This is why it is still a draft PR.

Checklist

• [ x] I added or updated tests for new behavior
• [ x] I updated documentation if user-facing behavior changed

See the developer guide for full coding standards.

GPU changes (expand if you modified src/simulation/)

• GPU results match CPU results
• Tested on NVIDIA GPU or AMD GPU

AI code reviews

Reviews are not retriggered automatically. To request a review, comment on the PR:

• @claude full review — Claude full review (also triggers on PR open/reopen/ready)
• Or add label claude-full-review — Claude full review via label

[github-actions]

Claude Code Review

Head SHA: 69cf365

Files changed:

• 16
• docs/documentation/case.md
• examples/2D_synthetic_turbulence/case.py
• src/common/m_constants.fpp
• src/common/m_global_parameters_common.fpp
• src/common/m_helper.fpp
• src/pre_process/m_perturbation.fpp
• src/simulation/m_body_forces.fpp
• src/simulation/m_checker.fpp
• src/simulation/m_global_parameters.fpp
• src/simulation/m_ibm.fpp
• src/simulation/m_mpi_proxy.fpp

Findings:

• src/simulation/m_checker.fpp:143-153 — s_check_inputs_synthetic_turbulence only checks num_turbulent_sources <= 0 and never checks it against num_turb_sources_max (=10, src/common/m_constants.fpp:119). turb_pos/synth_L are declared dimension(num_turb_sources_max, 3) (src/simulation/m_global_parameters.fpp:80). If a case sets num_turbulent_sources = 11, the checker's own loop do i = 1, num_turbulent_sources (line 145) reads turb_pos(i, d)/synth_L(i, d) out of bounds before the run even starts. The same gap exists for synth_n_shells vs num_synth_shells_max (=50): nothing prevents synth_n_shells from exceeding 50, and src/simulation/m_body_forces.fpp (num_synthetic_wave_numbers = sum(synth_n_waves_per_shell(1:synth_n_shells)) and the following do s = 1, synth_n_shells loop) will then index synth_k_shell/synth_amp_shell/synth_n_waves_per_shell past their declared size. This is a silent out-of-bounds array access (undetected on compilers without bounds-checking) rather than a clean error.

• src/simulation/m_body_forces.fpp:28,64,397 — num_synthetic_wave_numbers (module variable, no default initializer at line 28) is only assigned inside s_initialize_body_forces_module after the if (synth_n_shells <= 0) return guard at line 64. synth_n_shells defaults to dflt_int (src/simulation/m_global_parameters.fpp:484) and is not validated by the new checker, so a case with synthetic_turbulence = .true. and synth_n_shells left unset (or ≤ 0) hits this early return without ever assigning num_synthetic_wave_numbers. s_finalize_body_forces_module then reads this uninitialized value at line 397 (if (num_synthetic_wave_numbers > 0) then ... @:DEALLOCATE(...)), which can attempt to deallocate arrays that were never allocated — undefined behavior that will vary by compiler/build (debug vs release, gfortran vs nvfortran/Cray/ifx).

[sbryngelson]

Merged current master into this branch to resolve the conflict from #1550 (named parameter values / params codegen). One conflict in m_constants.fpp — your num_synth_shells_max/num_turb_sources_max constants kept alongside the new generated-constants include. No changes to your physics code.

Two pre-existing items will still fail CI, flagging so they're not a surprise:

1. The 5 new parameters (synthetic_turbulence, synth_seed, synth_n_shells, num_turbulent_sources, synth_U_inf) are registered but not documented in docs/documentation/case.md — precheck step 6 fails on that.
2. ./mfc.sh format wants changes in ~12 .fpp files on this branch (pre-dates the merge).

Both are quick fixes on your side. Also FYI: enumerated case parameters now accept named values ("riemann_solver": "hllc") — integer codes still work, no action needed.

[sbryngelson]

I've pushed a merge of master into this branch (5e64c8f, a true merge commit — your history is untouched) resolving the conflicts from the parameter-pipeline refactors (#1550–#1556). What was relocated and why:

src/simulation/m_global_parameters.fpp — your five synthetic-turbulence namelist scalars (synthetic_turbulence, synth_seed, synth_n_shells, num_turbulent_sources, synth_U_inf) are now declared by the auto-generated generated_decls.fpp (driven by your existing definitions.py registrations — those merged cleanly), so the manual declarations were removed to avoid duplicate-declaration errors. Their GPU_DECLARE moved to src/common/m_global_parameters_common.fpp (sim-only block): declare directives must live in the declaring module — Cray ftn rejects declare-target on use-associated names. The namelist arrays (synth_n_waves_per_shell, synth_k_shell, synth_amp_shell, turb_pos, synth_L) keep their manual declarations and GPU declares in the sim module (they're registered as indexed variants, which the generator deliberately leaves to the manual residue). Your defaults stay where they were; fft_wrt's default now comes from s_assign_common_defaults.

src/simulation/m_mpi_proxy.fpp — the big manual broadcast loops were replaced by the generated generated_bcast.fpp include, which now emits the broadcasts for all your scalars (and the nv_uvm_* ones) automatically. Only the five synth array broadcasts remain, in the manual-residue section with a comment explaining why.

I verified the generated output directly: all five scalars get declarations + broadcasts, all five arrays appear in the generated namelist, and the scalars are correctly absent from the manual residue. ./mfc.sh format/precheck, toolchain pytest (342 passed), and a full CPU build of all three targets pass on the merged tree.

For future parameters on this branch: registration in toolchain/mfc/params/definitions.py handles declaration/namelist/broadcast; only defaults (and GPU declares, in m_global_parameters_common) stay manual — see docs/documentation/contributing.md. Nothing left on your side; CI should run on the updated head.

[hyeoksu-lee]

Hi @danieljvickers @sbryngelson, you may already know, but I implemented a synthetic turbulence generator in pre_process/m_perturbation, specifically in s_perturb_mixlayer and s_generate_random_perturbation. This code comprises of two parts: 1) generates homogeneous turbulence based on a given spectrum, 2) transform it into anisotropic turbulence for mixing layer configuration. Part 1 should be more or less similar to what this PR introduces. This PR does similar tasks, but for different purpose and using different method, if I understood correctly.

One key difference would be my implementation runs only in the pre_process, because I targeted a temporal mixing layer, not a spatial one, while this PR is for a spatial configuration so the turbulence should be seeded at the inflow boundary over the simulation time.

Given that, I believe it would be better to put them together in a unified framework, if possible, in terms of DIY principles or at least for maintainability. Let me know if there's anything I can help you with.

[sbryngelson]

it certainly is a mess right now, and having a turbulence generator that's specific to a configuration (mixing layer) is niche

[sbryngelson]

@hyeoksu-lee - folded your generator in as the shared basis here.

• Moved s_prng, modmul, f_unit_vector from m_perturbation to src/common/m_helper (with $:GPU_ROUTINE); this PR's forcing now calls them instead of intrinsic random_number. s_perturb_mixlayer uses them unchanged via its existing use m_helper - mixlayer golden unaffected.
• Makes the forcing compiler-reproducible, so the 2D_synthetic_turbulence example is back in the CI golden suite (was skipped under intrinsic RNG).
• Used the double-cross-product polarization (sig = khat x (xi x khat)) for the 3-D solenoidal direction, fixing the non-divergence-free degenerate branch.

Left the call sites separate - yours is a pre_process IC with Reynolds-stress anisotropy, this is per-step inflow forcing. The spectrum/anisotropy layer could be unified later if worth it.

[sbryngelson]

Pushed two things to this branch and ran a validation pass on the synthetic-turbulence forcing (all on real MFC output, batch runs on hpcfund CPU).

Housekeeping pushed to the branch

• Merge conflict with master resolved (815f1e2f… now folded in) — the only conflict was the casesToSkip list; kept both 2D_synthetic_turbulence and master's new entries.
• New golden regression test (69cf3656): 3D -> synthetic_turbulence (CB853530). The feature had no automated coverage — the 2D_synthetic_turbulence example is in casesToSkip because the moving-airfoil-IB run is FP-marginal across compilers. This new test is a 24³ triply-periodic uniform single-fluid box driven by the deterministic 2-shell solenoidal forcing: IB-free and non-chaotic, so it holds across the 4 CI compilers at tolerance while exercising the full path (s_initialize_body_forces_module incl. the 3-D solenoidal construction, and s_compute_synthetic_forces_rhs). Generated and verified passing on CPU.

The two earlier Claude review findings are already fixed

Both were real at the reviewed SHA (dafeef62) but are resolved in 3d7c11d6 ("deterministic RNG + solenoidal polarization + mixture density"):

• Multi-fluid density — now sums the full cont%beg:cont%end range (mixture density), not just α₁ρ₁. ✅
• 3-D solenoidal direction — the buggy degenerate branch is gone; replaced by the double cross product σ = k̂ × (ξ × k̂), provably k·σ = 0. ✅

Validation results

1. Generator is exact — instrumented the mode generation: all 24 modes land on their prescribed shells, |σ| = 1, and max|k·σ| = 8×10⁻¹⁷ (machine-zero solenoidal — directly confirms the 3-D fix).
2. Assembled field — FFT of the injected velocity peaks at the three shells; the field is solenoidal (its discrete divergence matches a known-perfectly-solenoidal reference field).
3. Developed turbulence (128³, forced to t=3): Tu = u'/U∞ grows smoothly to ~0.5% and levels off; the developed energy spectrum peaks on the three forced shells with a cascade spreading to neighbors — i.e. the energy-distribution check you flagged as the remaining task.

One thing for the writeup

In the developed run, isotropy degrades (0.74 → 0.35; u'_x ends ~2.5× u'_y,z). With only ~24 modes plus the mean-flow x-advection the field isn't statistically isotropic — not a code issue, but the quantitative writeup would want more Fourier modes and/or ensemble/time averaging to demonstrate isotropy and a clean inertial range.

Happy to help wire up a larger isotropy/intensity study if useful.

[codecov]

Codecov Report

❌ Patch coverage is 64.91228% with 60 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.09%. Comparing base (102def2) to head (bfef557).
⚠️ Report is 3 commits behind head on master.

Files with missing lines	Patch %	Lines
src/simulation/m_body_forces.fpp	54.70%	29 Missing and 24 partials ⚠️
src/simulation/m_checker.fpp	50.00%	4 Missing and 1 partial ⚠️
src/common/m_helper.fpp	93.75%	1 Missing ⚠️
src/simulation/m_time_steppers.fpp	91.66%	0 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1548      +/-   ##
==========================================
+ Coverage   60.88%   61.09%   +0.20%     
==========================================
  Files          83       83              
  Lines       19836    20078     +242     
  Branches     2953     3002      +49     
==========================================
+ Hits        12077    12266     +189     
- Misses       5767     5783      +16     
- Partials     1992     2029      +37     

☔ 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.