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