fix: make cyl_to_geoflux consistent with cached forward transform

[krystophny]

User description

Summary

• Add invert_cached_surface() which uses Newton-Raphson iteration to find (s, θ) such that the cache interpolation maps back to a given (R, Z) point
• Modify cyl_to_geoflux and cyl_to_geoflux_internal to use cache-consistent inversion when cache is built
• Implement geoflux_from_cyl in the OOP interface to use the module's cyl_to_geoflux

Problem

The cyl_to_geoflux function was using direct psi evaluation while geoflux_to_cyl uses bilinear interpolation in the (s, θ) → (R, Z) cache. This caused roundtrip inconsistency because bilinear interpolation of (R, Z) between flux surface points gives points that are NOT on any flux surface (due to surface curvature).

Solution

When the cache is built, the inverse transformation now uses Newton-Raphson iteration to find the exact (s, θ) that maps through the cache interpolation to the given (R, Z) point. This makes the forward and inverse transformations mathematically consistent.

Test plan

• All existing tests pass (74/74)
• geoflux tests pass

---

PR Type

Bug fix, Enhancement

---

Description

• Add Newton-Raphson iteration for cache-consistent inverse transformation

• Implement invert_cached_surface() to find (s, θ) from (R, Z) coordinates

• Modify cyl_to_geoflux and cyl_to_geoflux_internal to use cache inversion

• Implement geoflux_from_cyl in OOP interface using module function

---

Diagram Walkthrough

flowchart LR
  A["(R, Z) input"] --> B["Cache built?"]
  B -->|Yes| C["invert_cached_surface<br/>Newton-Raphson iteration"]
  B -->|No| D["Direct psi evaluation"]
  C --> E["(s, θ) output"]
  D --> E
  F["geoflux_from_cyl<br/>OOP interface"] --> G["cyl_to_geoflux<br/>module function"]
  G --> E

Loading

File Walkthrough

	Relevant files
Bug fix	geoflux_coordinates.f90 Add cache-consistent inverse transformation with Newton-Raphson src/coordinates/geoflux_coordinates.f90 Add invert_cached_surface() subroutine using Newton-Raphson iteration with finite differences Modify cyl_to_geoflux() to use cache-consistent inversion when cache is built Modify cyl_to_geoflux_internal() to use cache-consistent inversion when cache is built Fallback to direct psi evaluation when cache is not available +61/-4
Enhancement	libneo_coordinates_geoflux.f90 Implement geoflux_from_cyl OOP interface method                    src/coordinates/libneo_coordinates_geoflux.f90 Import cyl_to_geoflux from geoflux_coordinates module Implement geoflux_from_cyl() subroutine using the module's cyl_to_geoflux function Replace error stop with functional implementation +7/-2

---

[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
⚪	🎫 No ticket provided Create ticket/issue
Codebase Duplication Compliance
⚪	Codebase context is not defined Follow the guide to enable codebase context checks.
Custom Compliance
🟢	Generic: Comprehensive Audit Trails Objective: To create a detailed and reliable record of critical system actions for security analysis and compliance. Status: Passed 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: Passed Learn more about managing compliance generic rules or creating your own custom rules
	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: Robust Error Handling and Edge Case Management Objective: Ensure comprehensive error handling that provides meaningful context and graceful degradation Status: Silent NR failure: invert_cached_surface can exit early on singular Jacobian or non-convergence and still returns potentially inaccurate (s_out, theta_out) without any status/diagnostic, making failures silent and hard to detect. Referred Code do iter = 1, max_iter call interpolate_cached_surface(s_cur, theta_cur, R_cur, Z_cur) err_R = R_cur - R_target err_Z = Z_cur - Z_target if (max(abs(err_R), abs(err_Z)) < tol_inv) exit call interpolate_cached_surface(clamp01(s_cur + h), theta_cur, R_sp, Z_sp) call interpolate_cached_surface(clamp01(s_cur - h), theta_cur, R_sm, Z_sm) dR_ds = (R_sp - R_sm) / (2.0_dp * h) dZ_ds = (Z_sp - Z_sm) / (2.0_dp * h) call interpolate_cached_surface(s_cur, theta_cur + h, R_tp, Z_tp) call interpolate_cached_surface(s_cur, theta_cur - h, R_tm, Z_tm) dR_dt = (R_tp - R_tm) / (2.0_dp * h) dZ_dt = (Z_tp - Z_tm) / (2.0_dp * h) det = dR_ds * dZ_dt - dR_dt * dZ_ds if (abs(det) < 1.0d-14) exit ... (clipped 10 lines) 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: Missing input bounds: invert_cached_surface accepts R_target/Z_target without local range validation (unlike cyl_to_geoflux which clamps), so callers such as cyl_to_geoflux_internal may pass out-of-domain values that could yield undefined interpolation behavior depending on interpolate_cached_surface. Referred Code if (ctx%cache_built) then call invert_cached_surface(xfrom(1), xfrom(3), s_val, theta_val) else s_val = s_from_position(xfrom(1), xfrom(3)) theta_val = atan2(xfrom(3) - ctx%Z_axis, xfrom(1) - ctx%R_axis) end if xto(1) = s_val xto(2) = theta_val xto(3) = xfrom(2) end subroutine cyl_to_geoflux_internal subroutine invert_cached_surface(R_target, Z_target, s_out, theta_out) real(dp), intent(in) :: R_target, Z_target real(dp), intent(out) :: s_out, theta_out real(dp) :: s_cur, theta_cur, R_cur, Z_cur real(dp) :: dR_ds, dR_dt, dZ_ds, dZ_dt real(dp) :: err_R, err_Z, det, ds, dt real(dp), parameter :: h = 1.0d-6 real(dp), parameter :: tol_inv = 1.0d-10 integer, parameter :: max_iter = 20 ... (clipped 36 lines) Learn more about managing compliance generic rules or creating your own custom rules

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	Handle solver non-convergence to prevent silent failures Add a check after the Newton-Raphson solver loop in invert_cached_surface to ensure convergence was achieved. If not, stop execution with an error to prevent silent failures and the propagation of inaccurate results. src/coordinates/geoflux_coordinates.f90 [816-842] do iter = 1, max_iter call interpolate_cached_surface(s_cur, theta_cur, R_cur, Z_cur) err_R = R_cur - R_target err_Z = Z_cur - Z_target if (max(abs(err_R), abs(err_Z)) < tol_inv) exit call interpolate_cached_surface(clamp01(s_cur + h), theta_cur, R_sp, Z_sp) call interpolate_cached_surface(clamp01(s_cur - h), theta_cur, R_sm, Z_sm) dR_ds = (R_sp - R_sm) / (2.0_dp * h) dZ_ds = (Z_sp - Z_sm) / (2.0_dp * h) call interpolate_cached_surface(s_cur, theta_cur + h, R_tp, Z_tp) call interpolate_cached_surface(s_cur, theta_cur - h, R_tm, Z_tm) dR_dt = (R_tp - R_tm) / (2.0_dp * h) dZ_dt = (Z_tp - Z_tm) / (2.0_dp * h) det = dR_ds * dZ_dt - dR_dt * dZ_ds if (abs(det) < 1.0d-14) exit ds = (dZ_dt * err_R - dR_dt * err_Z) / det dt = (-dZ_ds * err_R + dR_ds * err_Z) / det s_cur = clamp01(s_cur - ds) theta_cur = theta_cur - dt end do +if (iter > max_iter .or. abs(det) < 1.0d-14) then + error stop "invert_cached_surface: Failed to converge" +end if + Apply / Chat Suggestion importance[1-10]: 8 __ Why: This suggestion correctly identifies a critical silent failure mode where the solver may return inaccurate results without an error. Adding a convergence check and an error stop is a robust way to prevent the propagation of incorrect data in a physics simulation.	Medium
More