Add radial plasma profiles for collision calculations

[krystophny]

User description

Summary

Add support for radially-varying plasma profiles (density, temperature) enabling realistic collisionality calculations for fusion reactor conditions.

Features

• New profiles module with two profile types:

  • power_series: f(s) = scale * sum(coef(n) * s^n) - VMEC-compatible
  • two_power: f(s) = scale * (1 - s^p1)^p2 - supports non-integer exponents

• Collision coefficients precomputed on uniform s grid and interpolated at runtime

• Backward compatible: existing simulations without &profiles_nml section work unchanged

Usage Example

&profiles_nml
  profile_type = "two_power"
  
  Te_scale  = 10010.0d0   ! 10.01 keV on axis
  Te_p1 = 1.0d0
  Te_p2 = 2.0d0           ! Te = Te_scale * (1 - s)^2
  
  ni1_scale = 1.822d21    ! m^-3 on axis
  ni1_p1 = 1.0d0
  ni1_p2 = 3.5d0          ! ni1 = ni1_scale * (1 - s)^3.5
/

Test plan

• All 37 tests pass
• Power series evaluation tests
• Two-power evaluation tests (including non-integer exponents)
• Flat profile coefficients match scalar parameters
• Peaked profiles show higher collision rates in core vs edge

Closes #202

---

PR Type

Enhancement

---

Description

• Add new profiles module supporting radially-varying plasma profiles

  • Power series: f(s) = scale * sum(coef(n) * s^n)
  • Two-power: f(s) = scale * (1 - s^p1)^p2 with non-integer exponents

• Refactor collision calculations to precompute coefficients on uniform grid

• Implement linear interpolation for runtime coefficient lookup at arbitrary positions

• Maintain backward compatibility with existing simulations without profiles

• Add comprehensive unit tests for profile evaluation and collision rates

---

Diagram Walkthrough

flowchart LR
  A["profiles.f90<br/>Profile definitions"] --> B["init_collision_profiles<br/>Precompute grid"]
  B --> C["efcolf_grid<br/>Collision coefficients"]
  C --> D["get_local_coeffs<br/>Linear interpolation"]
  D --> E["stost<br/>Collision step"]
  F["simple_main<br/>read_profiles_config"] --> A

Loading

File Walkthrough

	Relevant files
Enhancement	profiles.f90 New profiles module with power series and two-power support src/profiles.f90 New module implementing two profile types: power_series and two_power Namelist-based configuration parsing for temperature and density profiles Pure functions for profile evaluation supporting non-integer exponents Central get_plasma_params subroutine dispatching to appropriate profile type +117/-0  collis_alphas.f90 Refactor collision calculations with profile grid support src/collis_alphas.f90 Refactored collision coefficient computation into separate compute_collision_coeffs subroutine Added coleff_local accepting collision coefficients as arguments for flexibility New init_collision_profiles precomputes coefficients on 101-point uniform s-grid New get_local_coeffs performs linear interpolation at arbitrary s positions Modified stost to use interpolated coefficients when profiles are enabled Modernized code style with explicit intent declarations and modern Fortran syntax +293/-298 simple_main.f90 Integrate profile initialization into main simulation        src/simple_main.f90 Added read_profiles_config subroutine calling read_profiles from profiles module Integrated profile configuration reading into main initialization sequence Modified init_collisions to call init_collision_profiles when profiles enabled Updated collide subroutine signature to accept inout z parameter +27/-6
Tests	test_profiles.f90 Add unit tests for profile evaluation and collision rates test/tests/test_profiles.f90 New comprehensive test suite with 4 test subroutines covering profile functionality Tests power series evaluation with linear and quadratic polynomials Tests two-power evaluation including non-integer exponents Validates flat profiles produce identical coefficients to scalar parameters Verifies peaked profiles show higher collision rates in core vs edge +268/-0
Configuration changes	CMakeLists.txt Add profiles module to build system                                            src/CMakeLists.txt Added profiles.f90 to source list before collis_alphas.f90 for proper dependency ordering +1/-0      CMakeLists.txt Add profile tests to test suite                                                    test/tests/CMakeLists.txt Added new test executable test_profiles.x linking against simple library Registered test with CMake test system with unit test label +5/-0

---

[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