Update initialization for turbulent mixing layer and add validation results

[hyeoksu-lee]

User description

Description

This PR includes an update on initialization method for turbulent mixing layer. Also, Matlab scripts for computing turbulence statistics are added to examples/3D_turb_mixing/turbulence_stat as well as validation results.

Type of change

Please delete options that are not relevant.

• Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Something else

Scope

• This PR comprises a set of related changes with a common goal

If you cannot check the above box, please split your PR into multiple PRs that each have a common goal.

How Has This Been Tested?

Ran examples/3D_turb_mixing on Carpenter

Streamwise velocity field at final time

Vorticity magnitude at final time

Evolution of Reynolds stress

Evolution of TKE budget

TKE budget averaged over self-similar period

Momentum thickness growth rate (see run_turbulence.m for references in legend)

Checklist

• I have added comments for the new code
• I added Doxygen docstrings to the new code
• I have made corresponding changes to the documentation (docs/)
• I have added regression tests to the test suite so that people can verify in the future that the feature is behaving as expected
• I have added example cases in examples/ that demonstrate my new feature performing as expected.
  They run to completion and demonstrate "interesting physics"
• I ran ./mfc.sh format before committing my code
• New and existing tests pass locally with my changes, including with GPU capability enabled (both NVIDIA hardware with NVHPC compilers and AMD hardware with CRAY compilers) and disabled
• This PR does not introduce any repeated code (it follows the DRY principle)
• I cannot think of a way to condense this code and reduce any introduced additional line count

---

PR Type

Enhancement, Documentation, Tests

---

Description

• Replace linear stability analysis with synthetic turbulence generation method for mixing layer initialization using Guo et al. (2023) approach
• Add comprehensive MATLAB scripts for turbulence statistics computation including Reynolds stress, TKE budget analysis, and validation against reference data
• Update 3D turbulent mixing layer example with higher Reynolds number (320), increased domain size and grid resolution
• Add new configuration parameters mixlayer_perturb_nk and mixlayer_perturb_k0 for synthetic turbulence control
• Remove deprecated mixlayer_domain parameter and associated complex validation requirements
• Implement pseudo-random number generator with supporting constants and helper functions
• Add documentation for new mixing layer perturbation method and turbulence analysis tools
• Update test cases and parameter dictionaries to support new synthetic turbulence approach

---

Changes walkthrough 📝

	Relevant files
Enhancement	6 files m_perturbation.fpp Replace linear stability analysis with synthetic turbulence generation src/pre_process/m_perturbation.fpp • Removed complex linear stability analysis code for mixing layer instability waves • Replaced s_superposition_instability_wave with new s_perturb_mixlayer function • Added spectrum-based synthetic turbulence generation method using Guo et al. (2023) approach • Implemented pseudo-random number generator and helper functions for perturbation generation +146/-527 m_initial_condition.fpp Update function call for mixing layer perturbation              src/pre_process/m_initial_condition.fpp • Changed function call from s_superposition_instability_wave to s_perturb_mixlayer +1/-1      m_constants.fpp Add constants for pseudo-random number generator                  src/common/m_constants.fpp • Added constants for pseudo-random number generator implementation • Defined modulus, multiplier, increment, amplifier, and decimal_trim parameters +7/-0      case.py Update 3D turbulent mixing layer example configuration      examples/3D_turb_mixing/case.py • Updated Reynolds number from 50 to 320 • Increased domain size and grid resolution significantly • Added grid stretching in y-direction with new parameters • Modified time stepping and output settings +21/-20  run_turbulence.m Add MATLAB script for turbulence statistics analysis          examples/3D_turb_mixing/turbulence_stat/run_turbulence.m • Added comprehensive MATLAB script for turbulence statistics computation • Implements Reynolds stress and TKE budget calculations • Includes plotting functions and comparison with reference data • Provides mixing layer thickness and growth rate analysis +601/-0  average_tke_over_self_similar.m Add TKE budget averaging and validation script                      examples/3D_turb_mixing/turbulence_stat/average_tke_over_self_similar.m • Added MATLAB script to compute time-averaged TKE budget over self-similar period • Implements interpolation and normalization of TKE transport, production, and dissipation terms • Includes plotting functionality with comparison to reference data from multiple studies • Contains helper functions for derivative computation and TKE budget visualization +118/-0
Configuration changes	7 files m_mpi_proxy.fpp Update MPI broadcasts for new mixing layer parameters        src/pre_process/m_mpi_proxy.fpp • Added MPI broadcast for new parameter mixlayer_perturb_nk • Added MPI broadcast for new parameter mixlayer_perturb_k0 • Removed MPI broadcast for deprecated parameter mixlayer_domain +4/-4      m_global_parameters.fpp Add new mixing layer perturbation parameters                          src/pre_process/m_global_parameters.fpp • Added new parameters mixlayer_perturb_nk and mixlayer_perturb_k0 for synthetic turbulence • Removed deprecated parameter mixlayer_domain • Set default values for new parameters (nk=100, k0=0.4446) +6/-2      m_checker.fpp Simplify mixing layer perturbation validation                        src/pre_process/m_checker.fpp • Simplified validation for mixlayer_perturb to only require p > 0 • Removed complex validation requirements for model equations, fluids, and boundary conditions +1/-5      m_start_up.fpp Update namelist for new mixing layer parameters                    src/pre_process/m_start_up.fpp • Updated namelist to include new parameters mixlayer_perturb_nk and mixlayer_perturb_k0 • Removed deprecated parameter mixlayer_domain from namelist +2/-2      case_dicts.py Update parameter dictionary for new mixing layer options  toolchain/mfc/run/case_dicts.py • Added new parameter types for mixlayer_perturb_nk and mixlayer_perturb_k0 • Removed deprecated parameter mixlayer_domain +2/-0      set_user_inputs.m Add user input configuration for turbulence statistics      examples/3D_turb_mixing/turbulence_stat/set_user_inputs.m • Added MATLAB function to configure simulation parameters for turbulent mixing layer analysis • Defines fluid properties, Reynolds number, grid configuration, and time stepping parameters • Sets up domain dimensions and stretched grid parameters for y-direction • Configures system indices for multi-fluid and multi-dimensional simulations +57/-0    .typos.toml Allow TKE acronym in typo checker                                                .typos.toml • Added TKE to the list of accepted words in typo checking configuration +1/-0
Tests	1 files cases.py Update test case for new mixing layer perturbation method toolchain/mfc/test/cases.py • Replaced complex alter_instability_wave function with simplified alter_mixlayer_perturb • Updated test case configuration for 3D mixing layer with new parameters • Simplified test case setup removing bubble-related configurations +33/-53
Miscellaneous	1 files golden-metadata.txt Update test metadata for new build environment                      tests/D13FDB23/golden-metadata.txt • Updated test metadata with new build configuration and environment • Changed from Delta cluster to M4Pro local machine • Updated compiler versions and CMake configuration +44/-68
Documentation	3 files case.md Update documentation for new mixing layer perturbation method docs/documentation/case.md • Updated documentation for mixing layer parameters • Removed mixlayer_domain parameter description • Added description of new synthetic turbulence generation method • Updated requirements for mixlayer_perturb option +3/-6      README.md Add README for turbulence statistics analysis tools            examples/3D_turb_mixing/turbulence_stat/README.md • Added README file describing turbulence statistics analysis tools • Provides references to relevant literature • Documents file structure and usage instructions +23/-0    references.md Add turbulence-related academic references                              docs/documentation/references.md • Added reference to Guo et al. (2023) for turbulence generation methods • Added reference to Michalke (1964) for hyperbolic-tangent velocity profile instability +4/-0
Additional files	9 files case.py +0/-2      reference.mat [link]    golden-metadata.txt +0/-174  golden.txt +0/-16    golden-metadata.txt +0/-174  golden.txt +0/-10    golden-metadata.txt +0/-174  golden.txt +0/-18    golden.txt +12/-12

---

Need help?

Type /help how to ... in the comments thread for any questions about Qodo Merge usage.

Check out the documentation for more information.

[codecov]

Codecov Report

Attention: Patch coverage is 76.05634% with 17 lines in your changes missing coverage. Please review.

Project coverage is 44.00%. Comparing base (4864d36) to head (c937bd9).
Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
src/pre_process/m_perturbation.fpp	75.00%	16 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #879      +/-   ##
==========================================
- Coverage   45.98%   44.00%   -1.98%     
==========================================
  Files          68       68              
  Lines       18629    18433     -196     
  Branches     2239     2228      -11     
==========================================
- Hits         8566     8112     -454     
- Misses       8711     9017     +306     
+ Partials     1352     1304      -48     

☔ View full report in Codecov by Sentry.
📢 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.

[hyeoksu-lee]

@sbryngelson Only frontier test fails. I guess something might be wrong with the pseudo-random number generator, but I am still trying to figure out what causes the failure. If you have any idea on what could be possibly wrong, that would be greatly appreciated!

[sbryngelson]

@hyeoksu-lee that is indeed weird. if I recall correctly, the Frontier CPU test passes? If that's correct, then it's very strange the GPU case would fail. I notice the failure is for a very small tolerance error. What does that particular test actually look like? Does it meaningfully change the state variables over the ~50 time steps (or whatever it is)? Perhaps it's just not a great test that is very sensitive to compiler optimizations or floating point roundoff.

[sbryngelson]

if you think the error in the Cray/Frontier case is due to the RNG, then you can print the random number it generates as a test and see it in the CI output (I think that should work, I believe we print the simulation output of failed tests in the CI output). also: Carpenter is a cray machine, so you can load cray modules there, compile MFC, and see if it passes tests!

[sbryngelson]

using Carpenter for this is probably the fastest way to debug. you'll have to change the modules from the ones ./mfc.sh load uses to the Cray ones. let me know if you can't figure this out, I have Carpenter access (i think)

[hyeoksu-lee]

@sbryngelson Thank you for the suggestions! I will try to verify with Carpenter's cray modules. Not sure if I can find all the required modules though. I'll keep you updated

[sbryngelson]

@sbryngelson Thank you for the suggestions! I will try to verify with Carpenter's cray modules. Not sure if I can find all the required modules though. I'll keep you updated

Let me know if that's the case and I'll try to help (put it in Slack).

[hyeoksu-lee]

@sbryngelson Using Carpenter Cray, I found that the PRNG was sensitive to machine precision errors due to the multiplication of a very large number followed by trimming. I’ve updated the PRNG to make it more robust against such precision issues. Test suite passed on both Carpenter GNU and Cray as well as my local macbook, so I believe the Frontier tests should now pass.

In the meantime, I will have to update the results of example/3D_turb_mixing as its initial condition has been changed and the statistics may change accordingly. I believe the statistics will remain close to the reference, but I want to make sure. This will take a day or two.

[sbryngelson]

@hyeoksu-lee ah, good find, very nice! Looks like it passed just fine. FYI the GT (Phoenix) runners are down or acting very strangely right now but I expect them to show fine results. I will mark this PR as a draft PR and you can mark it 'ready to review' once you check your final statistics.

[hyeoksu-lee]

@sbryngelson Statistics looks good. This PR is now ready for review.

[qodo-merge-pro]

PR Reviewer Guide 🔍

Here are some key observations to aid the review process:

⏱️ Estimated effort to review: 4 🔵🔵🔵🔵⚪

🧪 PR contains tests

🔒 No security concerns identified

⚡ Recommended focus areas for review Possible Issue The new perturbation method uses MPI broadcasts but only generates random numbers on rank 0. This could lead to synchronization issues or deadlocks if MPI communication fails, and the random number generation logic appears complex with potential for numerical instability. if (proc_rank == 0) then call s_generate_random_perturbation(khat, xi, phi, i, y_cc(r)) end if #ifdef MFC_MPI call MPI_BCAST(khat, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(xi, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(phi, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) #endif Code Quality The Cholesky decomposition implementation manually checks for small values and sets them to sgm_eps, but sgm_eps is not defined in the visible code. This could cause compilation errors or undefined behavior. if (abs(Lmat(1, 1)) < sgm_eps) Lmat(1, 1) = sgm_eps Lmat(2, 1) = Rij(2, 1)/Lmat(1, 1) Lmat(2, 2) = sqrt(Rij(2, 2) - Lmat(2, 1)**2._wp) if (abs(Lmat(2, 2)) < sgm_eps) Lmat(2, 2) = sgm_eps Lmat(3, 1) = Rij(3, 1)/Lmat(1, 1) Lmat(3, 2) = (Rij(3, 2) - Lmat(3, 1)*Lmat(2, 1))/Lmat(2, 2) Lmat(3, 3) = sqrt(Rij(3, 3) - Lmat(3, 1)**2._wp - Lmat(3, 2)**2._wp) Duplicate Code The derivative computation functions f_compute_derivative_1d and f_compute_derivative_3d contain very similar logic for boundary conditions and central differences. This could be refactored to reduce duplication. function dfunds = f_compute_derivative_1d(fun,s) dfunds = zeros(size(fun)); % initialize discrete derivative vector for j = 1:length(s) % for each index in s % if at bottom boundary, use 1-sided derivative dfunds(1) = (fun(2) - fun(1)) / (s(2) - s(1)); % if at top boundary, use 1-sided derivative dfunds(end) = (fun(end) - fun(end-1)) / (s(end) - s(end-1)); % otherwise, if in the interior of the domain, 2-sided derivative for i = 2:length(s)-1 dfunds(i) = (fun(i+1) - fun(i-1)) / (s(i+1) - s(i-1)); end end end % Compute the wall-normal derivative of a discretized function, fun(y) function dfunds = f_compute_derivative_3d(fun,s,dir) dfunds = zeros(size(fun)); % initialize discrete derivative vector if (dir == 1) % Forward difference at start dfunds(1,:,:) = (fun(2,:,:) - fun(1,:,:)) / (s(2) - s(1)); % Backward difference at end dfunds(end,:,:) = (fun(end,:,:) - fun(end-1,:,:)) / (s(end) - s(end-1)); % Central difference for interior for i = 2:length(s)-1 dfunds(i,:,:) = (fun(i+1,:,:) - fun(i-1,:,:)) / (s(i+1) - s(i-1)); end elseif (dir == 2) dfunds(:,1,:) = (fun(:,2,:) - fun(:,1,:)) / (s(2) - s(1)); dfunds(:,end,:) = (fun(:,end,:) - fun(:,end-1,:)) / (s(end) - s(end-1)); for i = 2:length(s)-1 dfunds(:,i,:) = (fun(:,i+1,:) - fun(:,i-1,:)) / (s(i+1) - s(i-1)); end elseif (dir == 3) dfunds(:,:,1) = (fun(:,:,2) - fun(:,:,1)) / (s(2) - s(1)); dfunds(:,:,end) = (fun(:,:,end) - fun(:,:,end-1)) / (s(end) - s(end-1)); for i = 2:length(s)-1 dfunds(:,:,i) = (fun(:,:,i+1) - fun(:,:,i-1)) / (s(i+1) - s(i-1)); end else disp("Wrong dir input"); end end

[Copilot]

Pull Request Overview

This PR replaces the previous linear-instability initialization of the turbulent mixing layer with a spectrum-based synthetic turbulence generator (Guo et al. 2023), updates all related configuration and broadcast parameters (mixlayer_perturb_nk, mixlayer_perturb_k0), and adds MATLAB tools (and documentation) for turbulence statistics and validation of the 3D mixing layer example.

• Replace s_superposition_instability_wave with s_perturb_mixlayer and remove legacy linear-stability code.
• Introduce new parameters mixlayer_perturb_nk and mixlayer_perturb_k0, propagate them through global parameters, MPI proxies, start-up, and namelists.
• Add MATLAB scripts (run_turbulence.m, average_tke_over_self_similar.m), example README, and reference documentation.

Reviewed Changes

Copilot reviewed 22 out of 31 changed files in this pull request and generated 1 comment.

Show a summary per file

File	Description
src/pre_process/m_perturbation.fpp	Remove linear-stability routines; add s_perturb_mixlayer subroutine and PRNG constants.
src/pre_process/m_initial_condition.fpp	Update initialization to call s_perturb_mixlayer instead.
src/common/m_constants.fpp	Define PRNG constants (modulus, multiplier, increment, etc.).
src/pre_process/m_global_parameters.fpp	Add mixlayer_perturb_nk and mixlayer_perturb_k0 defaults.
src/pre_process/m_mpi_proxy.fpp	Broadcast new mixing-layer parameters; remove deprecated mixlayer_domain.
src/pre_process/m_start_up.fpp	Include new perturbation parameters in namelist.
src/pre_process/m_checker.fpp	Simplify validation rules for mixlayer_perturb.
toolchain/mfc/run/case_dicts.py	Register new parameter types for mixlayer_perturb_nk and mixlayer_perturb_k0.
toolchain/mfc/test/cases.py	Replace legacy instability-wave test with alter_mixlayer_perturb.
examples/3D_turb_mixing/turbulence_stat/*	Add MATLAB scripts, README, and averaging tools for validation.
examples/3D_turb_mixing/case.py	Update 3D example to higher Re, grid stretching, and mixlayer config.
docs/documentation/case.md	Update docs to describe new synthetic turbulence parameters.
docs/documentation/references.md	Add Guo (2023) and Michalke (1964) references.
.typos.toml	Allow “TKE” in typo checker.

Comments suppressed due to low confidence (3)

examples/2D_mixing_artificial_Ma/case.py:110

• Enabling mixlayer_perturb for a 2D case (p == 0) will trigger the PROHIBIT check (p > 0) and cause startup to fail. Please disable or guard this setting in 2D examples.

            "mixlayer_vel_coef": 1.0,

src/pre_process/m_checker.fpp:191

• The new PROHIBIT rule only checks p > 0 but does not enforce mixlayer_vel_profile = .true.. To avoid misconfiguration, consider adding validation that the velocity profile flag is enabled when perturbations are requested.

        @:PROHIBIT(mixlayer_perturb .and. p == 0, "mixlayer_perturb requires p > 0")

toolchain/mfc/test/cases.py:634

• The new test generator only covers 3D (len(dimInfo[0]) == 3) cases. If 2D support is intended, add a branch to test mixlayer perturbations in 2D.

    def alter_mixlayer_perturb(dimInfo):

[hyeoksu-lee]

@sbryngelson I guess this is now ready for merging