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: dafeef6

Files changed:

• 10
• examples/2D_synthetic_turbulence/case.py
• src/common/m_constants.fpp
• src/simulation/m_body_forces.fpp
• src/simulation/m_global_parameters.fpp
• src/simulation/m_ibm.fpp
• src/simulation/m_mpi_proxy.fpp
• src/simulation/m_start_up.fpp
• src/simulation/m_time_steppers.fpp
• toolchain/mfc/params/definitions.py
• toolchain/mfc/params/descriptions.py

Findings

1. Wrong density used for force scaling in multi-fluid cases (correctness)

File: src/simulation/m_body_forces.fpp, s_compute_synthetic_forces_rhs

rho_local = q_prim_vf(eqn_idx%cont%beg)%sf(j, k, l)
f_scale = rho_local*(synth_U_inf/synth_L(turb_idx, 1))*G_norm

eqn_idx%cont%beg indexes the first partial density α₁ρ₁, not the mixture density. For num_fluids > 1 the forcing is scaled by a fraction of the total density. The PR's own example case has num_fluids=2 with alpha_rho(1)=0.8 and alpha_rho(2)=0.2, so the total density is 1.0 but rho_local would be 0.8.

The existing body-force path avoids this by using the precomputed rhoM array from s_compute_mixture_density, which sums all continuity components. The synthetic-force path should do the same or sum q_prim_vf(eqn_idx%cont%beg:eqn_idx%cont%end) directly.

---

2. Non-perpendicular solenoidal direction in 3-D degenerate branch (physics correctness)

File: src/simulation/m_body_forces.fpp, s_initialize_body_forces_module, degenerate 3-D branch

! k nearly parallel to z; cross with x-hat instead k_hat x x_hat = (0, -kz, k_y)/norm  (normalised)
synthetic_ex(m_global) = 0._wp
synthetic_ey(m_global) = -kz/max(sqrt(kz*kz + (k_mag*sin(rn1))**2), 1e-10_wp)
synthetic_ez(m_global) = (k_mag*sin(rn1))/max(sqrt(kz*kz + (k_mag*sin(rn1))**2), 1e-10_wp)

k_mag*sin(rn1) is used where r_xy*sin(rn1) is required. r_xy = sqrt(1-kz²) is the normalized xy-plane radius of the unit wave vector; k_mag is the physical wavenumber magnitude (e.g. 2π/0.6 ≈ 10.5). With k_mag substituted:

k · e = k_mag · sin(rn1) · kz · (k_mag − r_xy) / ‖e‖ ≠ 0

so e is NOT perpendicular to k, violating the divergence-free (solenoidal) constraint that is the stated purpose of this vector. The result is that each wave mode in the degenerate branch injects a compressible (acoustic) component in addition to the intended vortical forcing.

Fix: replace both occurrences of k_mag*sin(rn1) with r_xy*sin(rn1) in that branch:

synthetic_ey(m_global) = -kz/max(sqrt(kz*kz + (r_xy*sin(rn1))**2), 1e-10_wp)
synthetic_ez(m_global) = (r_xy*sin(rn1))/max(sqrt(kz*kz + (r_xy*sin(rn1))**2), 1e-10_wp)

This only affects 3-D runs where some wave vectors are nearly parallel to z (|kz| ≥ 0.9).

[wilfonba]

I think this module should be renamed if this is where you're going to add the synthetic turbulence code. The synthetic turbulence isn't really a body force since it's not applied everywhere in the domain. Though m_source_terms isn't a good name because there are other source terms elsewhere in the code. Maybe a standalone m_synthetic turbulence would be the best place to put this routine? There's nothing wrong with small modules. In many ways they're a good thing.

[wilfonba]

This line is relevant to my previous comment about code location

[wilfonba]

This line is relevant to my previous comment about code location

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