Add radial plasma profiles for collision calculations

[krystophny]

Risk tier

• T0: docs, comments, small build metadata
• T1: pure refactor, no behavior change intended
• T2: local numerical logic
• T3: physics, output behavior, coordinate convention
• T4: parallelism, GPU, dependency, CI, or security-sensitive build logic

Correctness contract

Intended behavior change

Add radially varying plasma profiles for collision calculations. The branch adds a profiles module, reads &profiles_nml, precomputes collision coefficients on a uniform s grid, and interpolates local coefficients during stochastic collision steps.

Behavior that must not change

Inputs without &profiles_nml keep the existing scalar collision behavior. Flat profiles must reproduce the scalar collision coefficients within the stated tolerance.

Coordinate / unit conventions

Profiles are functions of VMEC-like radial coordinate s. Temperatures, densities, and collision coefficients keep the existing units used by the collision module.

Numerical invariants

Flat radial profiles match scalar parameters. Peaked profiles produce larger core collision rates than edge collision rates. Non-integer two-power exponents evaluate without domain errors over the tested s interval.

Tests added

• unit: test_profiles.f90
• integration: focused test_profiles CTest target
• system: none
• golden record: unchanged in focused verification

Golden-record impact

• unchanged in focused verification
• changed

Failure modes considered

• Missing profile namelist must preserve old scalar behavior.
• Non-integer exponents at the plasma edge must not create invalid values.
• Interpolated coefficients must stay positive for physical profiles.
• This PR still needs the broader acceptance tests tracked in Define radial-collision physics acceptance tests #363 before merge: Maxwellian fixed point, pitch/energy moment relaxation, seeded stochastic reproducibility, and the relation to libneo collision/species ownership in Adopt libneo species and collision_freqs instead of own implementations #335.

Manual validation

Focused profile/collision coefficient tests passed after rebasing on current main.

Verification

$HOME/code/prompts/scripts/check-writing-slop.py src/profiles.f90 test/tests/test_profiles.f90
PASS: no writing-slop candidates at threshold medium

make test-nopy TEST=test_profiles VERBOSE=1
...
4:  Testing flat profile vs scalar coefficients...
4:    PASS: flat profile coefficients match scalar
4:  Testing peaked profile produces different collision rates...
4:    PASS: core has higher collision rate than edge
4:      core/edge ratio:   11.259823079577709
4:  All profile tests passed
1/1 Test #4: test_profiles ....................   Passed    0.09 sec
100% tests passed, 0 tests failed out of 1

[qodo-code-review]

PR Compliance Guide 🔍

Below is a summary of compliance checks for this PR:

Security Compliance
🟢	No security concerns identified No security vulnerabilities detected by AI analysis. Human verification advised for critical code.
Ticket Compliance
🟢	🎫 #202 🟢 Add a new profiles module to support radially varying plasma profiles (density and temperatures) as functions of s in [0,1]. Implement two analytical profile forms: power_series and two_power; allow non-integer exponents in two_power. Provide namelist configuration for profiles with on-axis scale factors and parameters; default to disabled/backward-compatible when absent or all zeros. Integrate profiles into collision calculations by precomputing collision coefficients on an s-grid and interpolating at runtime in stost using z(1) as s. Maintain backward compatibility by using existing scalar parameters when profiles are disabled. Ensure unit consistency (densities provided in SI m^-3 but collision routines use cm^-3; temperatures in eV). Initialize and read profiles during program startup before initializing collisions. Add tests to verify power_series and two_power evaluation, flat profile matches scalar behavior, and peaked profiles show higher collision rates in core than edge.
Codebase Duplication Compliance
⚪	Codebase context is not defined Follow the guide to enable codebase context checks.
Custom Compliance
🟢	Generic: Secure Error Handling Objective: To prevent the leakage of sensitive system information through error messages while providing sufficient detail for internal debugging. Status: Passed Learn more about managing compliance generic rules or creating your own custom rules
	Generic: Secure Logging Practices Objective: To ensure logs are useful for debugging and auditing without exposing sensitive information like PII, PHI, or cardholder data. Status: Passed Learn more about managing compliance generic rules or creating your own custom rules
⚪	Generic: Comprehensive Audit Trails Objective: To create a detailed and reliable record of critical system actions for security analysis and compliance. Status: No auditing: New runtime pathways (profile initialization, coefficient interpolation, collision updates) add critical simulation state changes without adding any structured audit logging of actions, parameters, user identity, or outcomes. Referred Code subroutine init_collision_profiles(am1, am2, Z1, Z2, ealpha, v0) use profiles, only: profiles_enabled, get_plasma_params real(wp), intent(in) :: am1, am2, Z1, Z2, ealpha, v0 real(wp) :: s, Te, Ti1, Ti2, ni1, ni2 real(wp) :: densi1_loc, densi2_loc, tempi1_loc, tempi2_loc, tempe_loc real(wp) :: v0_loc integer :: i if (.not. profiles_enabled) then use_profiles = .false. return end if use_profiles = .true. do i = 1, N_S_GRID s = dble(i - 1) / dble(N_S_GRID - 1) call get_plasma_params(s, Te, Ti1, Ti2, ni1, ni2) ... (clipped 32 lines) Learn more about managing compliance generic rules or creating your own custom rules
	Generic: Meaningful Naming and Self-Documenting Code Objective: Ensure all identifiers clearly express their purpose and intent, making code self-documenting Status: Vague names: Several short or cryptic identifiers (e.g., dp, dh, dpp, dhh, ur) persist in new/modified code and may hurt readability without local clarifying comments. Referred Code subroutine onseff(v, dp, dh, dpd) real(wp), intent(in) :: v real(wp), intent(out) :: dp, dh, dpd real(wp), parameter :: sqp = 1.7724538d0 real(wp), parameter :: cons = 0.75225278d0 real(wp) :: v2, v3, ex, er v2 = v**2 v3 = v2 * v if (v < 0.01d0) then dp = cons * (1.d0 - 0.6d0*v2) dh = cons * (1.d0 - 0.2d0*v2) dpd = 2.d0 * cons * (1.d0 - 1.2d0*v2) else if (v > 6.d0) then dp = 1.d0 / v3 dh = (1.d0 - 0.5d0/v2) / v dpd = -1.d0 / v3 else ex = exp(-v2) / sqp ... (clipped 6 lines) Learn more about managing compliance generic rules or creating your own custom rules
	Generic: Robust Error Handling and Edge Case Management Objective: Ensure comprehensive error handling that provides meaningful context and graceful degradation Status: Silent failure: read_profiles returns silently on open/read errors and only toggles profiles_enabled without surfacing actionable context or diagnostics, which could mask configuration issues. Referred Code subroutine read_profiles(config_file) character(len=*), intent(in) :: config_file integer :: ios open(unit=99, file=config_file, status='old', action='read', iostat=ios) if (ios /= 0) return read(99, nml=profiles_nml, iostat=ios) close(99) if (ios /= 0) then profiles_enabled = .false. return end if profiles_enabled = (profile_type /= "none") end subroutine read_profiles Learn more about managing compliance generic rules or creating your own custom rules
	Generic: Security-First Input Validation and Data Handling Objective: Ensure all data inputs are validated, sanitized, and handled securely to prevent vulnerabilities Status: Input validation: External configuration values (profile_type, scales, exponents) are read and used with minimal validation, which could lead to invalid states (e.g., negative densities/temperatures) despite some clamping later. Referred Code real(dp) :: Te_p1 = 1.0d0, Te_p2 = 0.0d0 real(dp) :: Ti1_p1 = 1.0d0, Ti1_p2 = 0.0d0 real(dp) :: Ti2_p1 = 1.0d0, Ti2_p2 = 0.0d0 real(dp) :: ni1_p1 = 1.0d0, ni1_p2 = 0.0d0 real(dp) :: ni2_p1 = 1.0d0, ni2_p2 = 0.0d0 logical :: profiles_enabled = .false. namelist /profiles_nml/ profile_type, & Te_scale, Ti1_scale, Ti2_scale, ni1_scale, ni2_scale, & Te_coef, Ti1_coef, Ti2_coef, ni1_coef, ni2_coef, & Te_p1, Te_p2, Ti1_p1, Ti1_p2, Ti2_p1, Ti2_p2, & ni1_p1, ni1_p2, ni2_p1, ni2_p2 contains Learn more about managing compliance generic rules or creating your own custom rules
Update

Compliance status legend

🟢 - Fully Compliant
🟡 - Partial Compliant
🔴 - Not Compliant
⚪ - Requires Further Human Verification
🏷️ - Compliance label

[qodo-code-review]

PR Code Suggestions ✨

Explore these optional code suggestions:

Category	Suggestion	Impact
Possible issue	Avoid division by zero in test Modify the test test_flat_profile_coefficients to prevent division by zero when calculating relative error by handling cases where the expected value is zero. test/tests/test_profiles.f90 [180-185] -if (abs(efcolf_profile(i) - efcolf_scalar(i)) / abs(efcolf_scalar(i)) > tol) then - print *, ' FAIL: efcolf mismatch for species', i - print *, ' scalar:', efcolf_scalar(i), 'profile:', efcolf_profile(i) - passed = .false. - nfail = nfail + 1 +if (abs(efcolf_scalar(i)) > tol) then + if (abs(efcolf_profile(i) - efcolf_scalar(i)) / abs(efcolf_scalar(i)) > tol) then + print *, ' FAIL: efcolf mismatch for species', i + print *, ' scalar:', efcolf_scalar(i), 'profile:', efcolf_profile(i) + passed = .false. + nfail = nfail + 1 + end if +else + if (abs(efcolf_profile(i) - efcolf_scalar(i)) > tol) then + print *, ' FAIL: efcolf mismatch for species', i + print *, ' scalar:', efcolf_scalar(i), 'profile:', efcolf_profile(i) + passed = .false. + nfail = nfail + 1 + end if end if Apply / Chat Suggestion importance[1-10]: 7 __ Why: The suggestion correctly identifies a potential division-by-zero error in the test suite if the expected scalar value is zero. The proposed fix of using a combined relative and absolute tolerance check makes the test more robust and prevents it from crashing.	Medium
	Prevent NaN from negative base exponentiation Add a check in eval_two_power to ensure s is non-negative before exponentiation to prevent potential NaN results. src/profiles.f90 [82-87] +if (s < 0.0d0) then + f = 0.0d0 + return +end if base = 1.0d0 - s**p1 if (base <= 0.0d0) then f = 0.0d0 else f = scale * base**p2 end if Apply / Chat Suggestion importance[1-10]: 5 __ Why: The suggestion correctly identifies a potential floating-point exception if s is negative and p1 is non-integer. Although s is a normalized coordinate and is expected to be non-negative, adding a defensive check improves the robustness of the pure function eval_two_power.	Low
Update