Add TF ripple perturbation to analytical tokamak field

[krystophny]

User description

Summary

Adds toroidal field (TF) coil ripple perturbation support to libneo's analytical tokamak equilibrium solver, following the ASCOT5 B_GS pattern.

Changes

Core Implementation

• TF Ripple Model: Added ripple parameters to analytical_circular_eq_t:
  • Nripple: Number of TF coils (0 = axisymmetric)
  • delta0: Ripple amplitude
  • alpha0, a0, z0: Ripple shape parameters
• Ripple Formula: δ(r,θ,φ) = δ₀(r/a₀)^α₀ exp(-θ²/2) cos(Nφ)
• Field Evaluation: New eval_bfield_ripple(R, φ, Z, ...) method with toroidal angle

Interface Extensions

• field_t Compatibility: Extended analytical_circular_eq_t to implement field_t interface
  • Enables use in Poincaré plots
  • Provides compute_bfield, compute_afield, compute_abfield methods
• Public Interface: Added compute_analytical_field_cylindrical() for external use

Testing

• Ripple Verification: Fortran test generates CSV data with/without ripple
• Automated Plotting: Python script creates 2D visualization via CTest
• field_t Test: Verifies polymorphic interface produces identical results
• Example Results: 9-coil configuration with δ₀=0.10 produces 12.65% peak-to-peak variation

Design Pattern

Following ASCOT5's B_GS.c approach:

• Single File: One type handles both axisymmetric and rippled fields
• Parameter-Controlled: Ripple is a parameter, not a separate type
• Dual Interface: Works as standalone and as field_t for Poincaré

Test Plan

• Ripple field generates expected periodic variation
• Automated plotting via CTest
• field_t interface matches native interface
• All existing tests pass

---

PR Type

Enhancement

---

Description

• Add TF ripple perturbation to analytical tokamak field

• Extend analytical_circular_eq_t to implement field_t interface

• Add comprehensive test suite with visualization

• Enable Poincaré plot compatibility

---

Diagram Walkthrough

flowchart LR
  A["analytical_circular_eq_t"] -- "extends" --> B["field_t interface"]
  A -- "adds" --> C["ripple parameters"]
  C -- "enables" --> D["eval_bfield_ripple()"]
  D -- "modulates" --> E["B_phi component"]
  F["test suite"] -- "generates" --> G["CSV data"]
  G -- "creates" --> H["2D visualization"]

Loading

File Walkthrough

	Relevant files
Enhancement	analytical_tokamak_field.f90 Add TF ripple and field_t interface                                            src/magfie/analytical_tokamak_field.f90 Add ripple parameters (Nripple, delta0, alpha0, a0, z0) to type Implement eval_bfield_ripple() method with toroidal angle support Extend type to inherit from field_t interface Add field_t interface methods (compute_bfield, compute_afield, compute_abfield) Add public compute_analytical_field_cylindrical() function +162/-6
Tests	test_analytical_circular.f90 Add field_t interface test                                                              test/source/test_analytical_circular.f90 Add test_field_t_interface() subroutine Verify field_t interface produces identical results to native interface Test polymorphic usage via field_t pointer +53/-0    test_ripple_field.f90 Add ripple field test program                                                        test/source/test_ripple_field.f90 New test program for TF ripple functionality Generate CSV data for ripple comparison (with/without ripple) Test 9-coil configuration with δ₀=0.10 ripple amplitude Output data for visualization pipeline +92/-0    plot_ripple_comparison.py Add ripple visualization script                                                    test/scripts/plot_ripple_comparison.py New Python script for ripple visualization Create 2D pcolormesh plots in (φ,Z) coordinates Show unperturbed, perturbed, and difference fields Calculate ripple statistics and periodicity verification +124/-0  CMakeLists.txt Add ripple test configuration                                                        test/CMakeLists.txt Add test_ripple_field.x executable target Add test_ripple_field_data test for data generation Add test_ripple_field_plot test for automated visualization Set up test dependencies between data generation and plotting +17/-0

---

[qodo-code-review]

PR Compliance Guide 🔍

Below is a summary of compliance checks for this PR:

Security Compliance
⚪	Unimplemented method abort Description: The method compute_afield is declared but intentionally unimplemented and aborts execution via error stop, which can cause unexpected termination if called through the field_t interface. analytical_tokamak_field.f90 [371-372] Referred Code error stop 'analytical_circular_eq_t: compute_afield not implemented' end subroutine compute_afield
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
⚪	No custom compliance provided Follow the guide to enable custom compliance check.
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	Ensure magnetic field is divergence-free To ensure the magnetic field remains divergence-free, add perturbations to the radial (B_R) and vertical (B_Z) components derived from a scalar potential, instead of only modifying the toroidal component (B_phi). src/magfie/analytical_tokamak_field.f90 [301-331] subroutine eval_bfield_ripple(self, R, phi, Z, B_R, B_Z, B_phi, B_mod) class(analytical_circular_eq_t), intent(in) :: self real(dp), intent(in) :: R, phi, Z real(dp), intent(out) :: B_R, B_Z, B_phi, B_mod real(dp) :: dpsi_dR, dpsi_dZ real(dp) :: F_psi real(dp) :: radius, theta, delta_ripple + real(dp) :: d_delta_dR, d_delta_dZ, d_radius_dR, d_radius_dZ, d_theta_dR, d_theta_dZ + real(dp) :: delta_B_R, delta_B_Z + real(dp), parameter :: tiny = 1.0e-12_dp if (.not. self%initialized) then error stop "analytical_circular_eq_t not initialized" end if call self%eval_psi_derivatives(R, Z, dpsi_dR, dpsi_dZ) B_R = -(1.0_dp / R) * dpsi_dZ B_Z = (1.0_dp / R) * dpsi_dR F_psi = self%B0 * self%R0 B_phi = F_psi / R if (self%Nripple > 0) then radius = sqrt((R - self%R0)**2 + (Z - self%z0)**2) theta = atan2(Z - self%z0, R - self%R0) delta_ripple = self%delta0 * exp(-0.5_dp * theta**2) & * (radius / self%a0)**self%alpha0 B_phi = B_phi * (1.0_dp + delta_ripple * cos(real(self%Nripple, dp) * phi)) + + ! Add divergence-free ripple components to B_R and B_Z + if (radius > tiny) then + d_radius_dR = (R - self%R0) / radius + d_radius_dZ = (Z - self%z0) / radius + d_theta_dR = -(Z - self%z0) / radius**2 + d_theta_dZ = (R - self%R0) / radius**2 + + d_delta_dR = delta_ripple * ((-theta * d_theta_dR) + (self%alpha0 / radius * d_radius_dR)) + d_delta_dZ = delta_ripple * ((-theta * d_theta_dZ) + (self%alpha0 / radius * d_radius_dZ)) + + delta_B_R = (self%B0 * self%R0 / real(self%Nripple, dp)) * d_delta_dR * sin(real(self%Nripple, dp) * phi) + delta_B_Z = (self%B0 * self%R0 / real(self%Nripple, dp)) * d_delta_dZ * sin(real(self%Nripple, dp) * phi) + + B_R = B_R + delta_B_R + B_Z = B_Z + delta_B_Z + end if end if B_mod = sqrt(B_R**2 + B_Z**2 + B_phi**2) end subroutine eval_bfield_ripple Apply / Chat Suggestion importance[1-10]: 9 __ Why: The suggestion correctly identifies that only perturbing B_phi violates the divergence-free condition for the magnetic field, which is a critical physical constraint, and proposes a valid correction.	High
High-level	Unify B-field evaluation methods Unify the new eval_bfield_ripple and existing eval_bfield methods into a single B-field calculation routine. The new method already handles the non-rippled case, making the old one redundant. Examples: src/magfie/analytical_tokamak_field.f90 [289-331] !> Evaluate magnetic field components with TF ripple !> !> Parameters: !> R - Major radius coordinate [m] !> phi - Toroidal angle [rad] !> Z - Vertical coordinate [m] !> !> Returns: !> B_R - Radial field component [T] !> B_Z - Vertical field component [T] ... (clipped 33 lines) Solution Walkthrough: Before: type :: analytical_circular_eq_t ... contains procedure :: eval_bfield => evaluate_bfield procedure :: eval_bfield_ripple end type subroutine evaluate_bfield(self, R, Z, B_R, B_Z, B_phi, B_mod) ! ... calculate axisymmetric B-field ... end subroutine subroutine eval_bfield_ripple(self, R, phi, Z, B_R, B_Z, B_phi, B_mod) ! ... calculate axisymmetric B-field (duplicated logic) ... if (self%Nripple > 0) then ! ... add ripple to B_phi ... end if end subroutine After: type :: analytical_circular_eq_t ... contains procedure :: eval_bfield => evaluate_bfield end type subroutine evaluate_bfield(self, R, Z, B_R, B_Z, B_phi, B_mod, phi) real(dp), intent(in), optional :: phi ! ... calculate axisymmetric B-field ... if (self%Nripple > 0 .and. present(phi)) then ! ... add ripple to B_phi ... end if B_mod = sqrt(B_R**2 + B_Z**2 + B_phi**2) end subroutine Suggestion importance[1-10]: 7 __ Why: The suggestion correctly identifies significant code duplication between eval_bfield and eval_bfield_ripple, and unifying them would improve maintainability and simplify the class interface.	Medium
General	Improve flux surface finding accuracy Replace the coarse grid search for finding a flux surface in the test with a more robust and accurate bisection method to improve test reliability. test/source/test_ripple_field.f90 [43-58] found_flux_surface = .false. -do i = 1, nR - R_flux = R_min + (R_max - R_min) * real(i - 1, dp) / real(nR - 1, dp) - psi_val = eq_no_ripple%eval_psi(R_flux, 0.0_dp) - if (abs(psi_val - psi_target) < 1.0e-3_dp * abs(psi_target) .and. R_flux > R0) then - found_flux_surface = .true. - exit - end if -end do +block + real(dp) :: R_low, R_high, R_mid, psi_mid, psi_low + real(dp), parameter :: tol = 1.0e-6_dp + integer :: k, max_iter + max_iter = 100 + R_low = R0 + R_high = R_max + psi_low = eq_no_ripple%eval_psi(R_low, 0.0_dp) + + do k = 1, max_iter + R_mid = 0.5_dp * (R_low + R_high) + psi_mid = eq_no_ripple%eval_psi(R_mid, 0.0_dp) + if (abs(R_high - R_low) < tol) then + found_flux_surface = .true. + R_flux = R_mid + exit + end if + if (sign(1.0_dp, psi_mid - psi_target) == sign(1.0_dp, psi_low - psi_target)) then + R_low = R_mid + psi_low = psi_mid + else + R_high = R_mid + end if + end do +end block if (.not. found_flux_surface) then R_flux = R0 + 0.5_dp - print *, "Warning: Could not find exact flux surface, using R =", R_flux + print *, "Warning: Bisection failed, using R =", R_flux else print *, "Found flux surface at R =", R_flux end if Apply / Chat Suggestion importance[1-10]: 7 __ Why: The suggestion improves the test's robustness by replacing a simple grid search with a more accurate bisection method for finding the flux surface, leading to more reliable test data.	Medium
	Avoid hardcoding parameters in plots Avoid hardcoding ripple parameters in the plot title. Instead, pass them from the Fortran test to the Python plotting script to generate the title dynamically, ensuring it matches the data. test/scripts/plot_ripple_comparison.py [85-91] +# At the top of main(), after parsing args: +# ... +# params_file = data_dir / 'ripple_params.txt' +# with open(params_file, 'r') as f: +# n_ripple = int(f.readline().split('=')[1]) +# delta0 = float(f.readline().split('=')[1]) +# ... + +# In the plotting section: im1 = axes[1].pcolormesh(phi_deg_rip, Z_grid_rip, Bphi_grid_rip, shading='auto', vmin=vmin, vmax=vmax, cmap='viridis') axes[1].set_xlabel('Toroidal angle φ [deg]', fontsize=12) axes[1].set_ylabel('Z [m]', fontsize=12) -axes[1].set_title('B_φ with 9-coil ripple (δ₀=0.10)', fontsize=14) +axes[1].set_title(f'B_φ with {n_ripple}-coil ripple (δ₀={delta0:.2f})', fontsize=14) axes[1].set_aspect('auto') plt.colorbar(im1, ax=axes[1], label='B_φ [T]') [To ensure code accuracy, apply this suggestion manually] Suggestion importance[1-10]: 5 __ Why: The suggestion correctly points out that hardcoded parameters in the plot title are brittle and proposes a good solution to make the test visualization more robust by dynamically generating the title.	Low
Update