Add TF ripple perturbation to analytical tokamak field#149
Conversation
PR Compliance Guide 🔍Below is a summary of compliance checks for this PR:
Compliance status legend🟢 - Fully Compliant🟡 - Partial Compliant 🔴 - Not Compliant ⚪ - Requires Further Human Verification 🏷️ - Compliance label |
||||||||||||||||||
krystophny
temporarily deployed
to
github-pages
GitHub Actions
Inactive
PR Code Suggestions ✨Explore these optional code suggestions:
|
|||||||||||||||||
krystophny
temporarily deployed
to
github-pages
GitHub Actions
Inactive
krystophny
temporarily deployed
to
github-pages
GitHub Actions
Inactive
Implements toroidal field ripple following ASCOT5's B_GS model with δ(r,θ,φ) = δ₀ (r/a₀)^α₀ exp(-θ²/2) cos(N_ripple φ) Changes: - Add ripple parameters to analytical_circular_eq_t type - Nripple: number of TF coils (0 = no ripple) - delta0: ripple amplitude - alpha0: radial dependency exponent - a0: minor radius reference - z0: midplane reference - New eval_bfield_ripple method accepts toroidal angle phi - Ripple modulates B_phi component while preserving poloidal field - Full test suite with automatic plotting via CTest Tests verify: - 9-coil ripple with δ₀=0.10 produces ~12.65% peak-to-peak variation - 40° periodicity (360°/9) visible in 2D (φ,Z) pcolormesh plots - Axisymmetric field recovered when Nripple=0 Test files: - test/source/test_ripple_field.f90: Fortran data generation - test/scripts/plot_ripple_comparison.py: Python visualization - Outputs: bphi_ripple_comparison.png showing unperturbed, perturbed, and difference fields
Following ASCOT5's B_GS pattern, the analytical equilibrium type now extends field_t directly instead of requiring a separate wrapper file. This allows the same type to be used both as a standalone field (via eval_bfield_ripple) and as a Poincare-compatible field_t. Added methods: - compute_bfield: field_t interface for magnetic field evaluation - compute_afield: field_t interface stub (not implemented) - compute_abfield: combined A and B evaluation - compute_analytical_field_cylindrical: public standalone function The Nripple parameter controls ripple: 0 = axisymmetric, >0 = rippled.
Tests that analytical_circular_eq_t can be used as a field_t (via pointer) and produces identical results to the direct eval_bfield_ripple interface. This verifies Poincare plot compatibility. The test creates a field_t pointer to an analytical equilibrium instance, evaluates the magnetic field at a test point using both the field_t interface and the native interface, and verifies they produce identical results within numerical tolerance (1e-12).
- Unit tests: analytical_circular, analytical_geoflux, ripple_field - Integration tests: analytical_geoflux_integration - All tests pass (3/3) - Documents field-agnostic geoflux achievement - Outlines next steps for SIMPLE system tests
krystophny
force-pushed
the
feature/tf-ripple-perturbation
branch
from
93e0cec to
41a628f
Compare
krystophny
temporarily deployed
to
github-pages
GitHub Actions
Inactive
eval_bfield_ripple previously perturbed only B_phi by B_phi*(1+delta0*cos(N*phi)), with B_R/B_Z left axisymmetric. This violates div(B)=0 by O(delta0) (verified up to 2.6% relative divergence via finite differences). Replace it with an exactly divergence-free construction: add the ripple as curl(A_Z_ripple e_Z) with A_Z_ripple = -(B0*R0/Nripple)*delta(R,Z)*cos(N*phi), giving B_R += (B0*R0/R) * delta(R,Z) * sin(N*phi) B_phi += (B0*R0/Nripple) * d(delta)/dR * cos(N*phi) B_Z unchanged which is solenoidal by construction (curl is always divergence free) for any envelope delta(R,Z), verified to ~1e-10 relative residual by finite differences over a grid of (R, phi, Z) points. compute_afield unconditionally did `error stop`, so the public compute_abfield (which calls it) always crashed. Implement it consistently with the new field: A = (0, psi/R, A_Z) with A_Z = -B0*R0*ln(R) - (B0*R0/Nripple)*delta(R,Z)*cos(N*phi), verified by finite-difference curl(A) against compute_bfield for both Nripple=0 and Nripple>0. Replace the CSV-dump/PNG-plot-only ripple "tests" with real assertions: N-fold periodicity (and non-periodicity at half that period), peak-to-peak B_phi modulation of the same order as delta0, linear scaling with delta0, and exact recovery of the axisymmetric field at Nripple=0. Delete TEST_SUMMARY.md, which documented tests/results that don't exist in the codebase.
![slopqueue[bot]](https://avatars.githubusercontent.com/u/48486978?s=60&v=4){kind=small}
There was a problem hiding this comment.
ORCHESTRATE_REVIEW_STATUS: approve
The PR adds TF ripple perturbation to the analytical tokamak field via a divergence-free curl(A) construction, extends the type to implement field_t for Poincaré compatibility, and adds comprehensive tests. The physics is correct: B = curl(A) is exact (verified by manual calculation and the test_compute_abfield_curl_consistency FD test), div(B) = 0 is verified numerically, and the coordinate convention x=(R,φ,Z), B=(B_R,B_φ,B_Z) matches circular_tokamak_field and the Poincaré consumer in src/poincare.f90:33-41. All CI checks pass (Debug, Release, docs, gate). The magfie library correctly links field_base transitively through neo. Findings below are all non-blocking.
Finding 1 — Unused pi parameter (compiler warning)
test/source/test_analytical_circular.f90:409: real(dp), parameter :: pi = ... is declared in test_divergence_free_ripple but never used. The Debug CI build emits Warning: Unused parameter 'pi' [-Wunused-parameter]. Fix: remove the line.
Finding 2 — ddelta_dR computed but unused in compute_afield
src/magfie/analytical_tokamak_field.f90:437: ripple_envelope returns both delta_ripple and ddelta_dR, but compute_afield only uses delta_ripple. The derivative is wasted computation. Not a bug; consider a delta_only helper if this matters for hot paths.
Finding 3 — test_ripple_amplitude_order_delta0 is very loose
test/source/test_ripple_field.f90:83-110: Checks that peak-to-peak B_φ ripple is within [0.3·δ₀, 3·δ₀] — a 10x window. The companion test_ripple_amplitude_scales_with_delta0 (exact 2× scaling, tol 1e-6) compensates. Non-blocking, but the loose test adds little beyond "doesn't crash."
Finding 4 — test_field_t_interface is tautological
test/source/test_analytical_circular.f90:471-510: Compares compute_bfield (which delegates to eval_bfield_ripple) against eval_bfield_ripple directly. Only verifies component mapping B=(B_R,B_φ,B_Z), not field values. The test_compute_abfield_curl_consistency test independently verifies field values via curl, so coverage is adequate.
Finding 5 — compute_analytical_field_cylindrical untested directly
src/magfie/analytical_tokamak_field.f90:398-404: Public one-line wrapper around eval_bfield_ripple. Tested indirectly. No external callers exist in the repo.
Finding 6 — No validation of edge-case ripple parameters
src/magfie/analytical_tokamak_field.f90:296-318: a0=0 → division by zero in ripple_envelope; alpha0 ≤ 1 → regularization at radius=0 is incorrect (comment says "for alpha0 > 1"). Defaults are safe (a0=R₀·ε, alpha0=2.0). Non-blocking since these require explicit user action.
Finding 7 — PR body describes files not in the diff
The PR body references plot_ripple_comparison.py, test_ripple_field_data, and test_ripple_field_plot, none of which appear in the actual 4-file diff. The auto-generated walkthrough is stale.
{"findings":[{"file":"test/source/test_analytical_circular.f90","line":409,"severity":"low","message":"Unused 'pi' parameter triggers -Wunused-parameter warning in Debug CI"},{"file":"src/magfie/analytical_tokamak_field.f90","line":437,"severity":"low","message":"ddelta_dR computed by ripple_envelope but unused in compute_afield — wasted computation"},{"file":"test/source/test_ripple_field.f90","line":83,"severity":"low","message":"test_ripple_amplitude_order_delta0 uses a 10x tolerance window (0.3x-3x delta0); very loose, compensated by scaling test"},{"file":"test/source/test_analytical_circular.f90","line":471,"severity":"low","message":"test_field_t_interface is tautological: compares compute_bfield with eval_bfield_ripple it delegates to; only tests component mapping"},{"file":"src/magfie/analytical_tokamak_field.f90","line":398,"severity":"low","message":"compute_analytical_field_cylindrical is public but not tested directly; one-line wrapper"},{"file":"src/magfie/analytical_tokamak_field.f90","line":296,"severity":"low","message":"No validation of edge-case parameters: a0=0 causes division by zero, alpha0<=1 makes regularization incorrect; defaults are safe"}]}| real(dp) :: div_B, tol | ||
| integer :: i, j, k | ||
| real(dp) :: r_min, r_max, z_max | ||
| real(dp), parameter :: pi = 4.0_dp * atan(1.0_dp) |
There was a problem hiding this comment.
low: Unused 'pi' parameter triggers -Wunused-parameter warning in Debug CI
| A(3) = -F_psi * log(R) | ||
|
|
||
| if (self%Nripple > 0) then | ||
| call ripple_envelope(self, R, Z, delta_ripple, ddelta_dR) |
There was a problem hiding this comment.
low: ddelta_dR computed by ripple_envelope but unused in compute_afield — wasted computation
| print *, "" | ||
| end subroutine test_ripple_periodicity | ||
|
|
||
| subroutine test_ripple_amplitude_order_delta0() |
There was a problem hiding this comment.
low: test_ripple_amplitude_order_delta0 uses a 10x tolerance window (0.3x-3x delta0); very loose, compensated by scaling test
| subroutine test_field_t_interface() | ||
| !> Test field_t interface (for Poincare compatibility) | ||
| use neo_field_base, only: field_t | ||
| type(analytical_circular_eq_t), target :: eq |
There was a problem hiding this comment.
low: test_field_t_interface is tautological: compares compute_bfield with eval_bfield_ripple it delegates to; only tests component mapping
| !> B_Z - Vertical field component [T] | ||
| !> B_phi - Toroidal field component (with ripple if enabled) [T] | ||
| !> B_mod - Total field magnitude [T] | ||
| subroutine compute_analytical_field_cylindrical(eq, R, phi, Z, B_R, B_Z, B_phi, B_mod) |
There was a problem hiding this comment.
low: compute_analytical_field_cylindrical is public but not tested directly; one-line wrapper
| !> ddelta_dR = d(delta)/dR at fixed Z, regularized to zero at radius=0 | ||
| !> where the (radius, theta) parametrization is singular but delta and | ||
| !> its gradient vanish in the limit for alpha0 > 1. | ||
| subroutine ripple_envelope(self, R, Z, delta_ripple, ddelta_dR) |
There was a problem hiding this comment.
low: No validation of edge-case parameters: a0=0 causes division by zero, alpha0<=1 makes regularization incorrect; defaults are safe
1 participant
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
analytical_circular_eq_t:Nripple: Number of TF coils (0 = axisymmetric)delta0: Ripple amplitudealpha0,a0,z0: Ripple shape parameterseval_bfield_ripple(R, φ, Z, ...)method with toroidal angleInterface Extensions
analytical_circular_eq_tto implementfield_tinterfacecompute_bfield,compute_afield,compute_abfieldmethodscompute_analytical_field_cylindrical()for external useTesting
Design Pattern
Following ASCOT5's B_GS.c approach:
Test Plan
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
File Walkthrough
analytical_tokamak_field.f90
Add TF ripple and field_t interfacesrc/magfie/analytical_tokamak_field.f90
compute_abfield)
test_analytical_circular.f90
Add field_t interface testtest/source/test_analytical_circular.f90
interface
test_ripple_field.f90
Add ripple field test programtest/source/test_ripple_field.f90
plot_ripple_comparison.py
Add ripple visualization scripttest/scripts/plot_ripple_comparison.py
CMakeLists.txt
Add ripple test configurationtest/CMakeLists.txt