docs: Add GEQDSK/geoflux implementation plan with 3D field support

CMakeLists.txt

@@ -46,6 +46,11 @@ option(ENABLE_GVEC "Enable GVEC field support (experimental)" OFF)
 option(ENABLE_COVERAGE "Enable code coverage analysis (Debug/Profile builds only)" OFF)
 option(SIMPLE_DETERMINISTIC_FP "Disable fast-math for reproducible floating-point" OFF)
 
+# Golden record regression tests require deterministic floating-point output.
+if (SIMPLE_TESTING)
+    set(SIMPLE_DETERMINISTIC_FP ON CACHE BOOL "Disable fast-math for reproducible floating-point" FORCE)
+endif()
+
 # Conditionally fetch GVEC if enabled
 if(ENABLE_GVEC)
     message(STATUS "GVEC support enabled - fetching GVEC library")
@@ -92,8 +97,53 @@ add_link_options(${NETCDF_FLIBS})
 find_package(BLAS REQUIRED)
 find_package(LAPACK REQUIRED)
 
+if (SIMPLE_TESTING)
+    # Golden record tests (and CI) must not silently pick up local checkouts via
+    # the CODE environment variable, otherwise the current build can differ from
+    # the reference build even without code changes.
+    if (DEFINED ENV{CODE} AND NOT "$ENV{CODE}" STREQUAL "")
+        message(STATUS "SIMPLE_TESTING: ignoring CODE=$ENV{CODE} for reproducible dependency resolution")
+        set(ENV{CODE} "")
+    endif()
+endif()
+
 find_or_fetch(libneo)
 
+# libneo enables -ffast-math in Release/Fast by default, which breaks strict
+# reproducibility. In deterministic mode, strip fast-math options from the libneo
+# targets that SIMPLE links against.
+if (SIMPLE_DETERMINISTIC_FP AND CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
+    function(_simple_strip_libneo_fast_math target_name)
+        if (NOT TARGET ${target_name})
+            return()
+        endif()
+
+        get_target_property(_opts ${target_name} COMPILE_OPTIONS)
+        if (NOT _opts)
+            target_compile_options(${target_name} PRIVATE -ffp-contract=off)
+            return()
+        endif()
+
+        set(_new_opts "")
+        foreach(_opt IN LISTS _opts)
+            if (_opt STREQUAL "-ffast-math")
+                continue()
+            endif()
+            if (_opt STREQUAL "-ffp-contract=fast")
+                continue()
+            endif()
+            list(APPEND _new_opts "${_opt}")
+        endforeach()
+
+        set_target_properties(${target_name} PROPERTIES COMPILE_OPTIONS "${_new_opts}")
+        target_compile_options(${target_name} PRIVATE -ffp-contract=off)
+    endfunction()
+
+    foreach(_libneo_target neo neo_field interpolate neo_polylag hdf5_tools CONTRIB magfie)
+        _simple_strip_libneo_fast_math(${_libneo_target})
+    endforeach()
+endif()
+
 # Fetch fortplot for diagnostic plotting
 include(FetchContent)
 FetchContent_Declare(

TODO.md

@@ -0,0 +1,136 @@
+# TODO: Meiss Canonical Coordinates on EQDSK/Geoflux — Full Implementation & Robust Testing
+
+_Last updated: 2025-09-30_
+
+This is an exhaustive execution guide. Every step includes file names, code snippets, commands, and validation criteria. The emphasis is on **precise, fast, non-tautological tests** and **system-level regressions**.
+
+---
+## 0. Baseline Safety Net (Tests Only)
+1. **Ensure clean trees** (SIMPLE + libneo):
+   ```bash
+   git status
+   cd ../libneo && git status
+   cd ../SIMPLE
+   ```
+
+2. **VMEC Meiss guard** (`test/tests/test_field_can_meiss_vmec.f90`)
+   - Implement as described previously **but add analytic invariants**:
+     - After `evaluate_meiss`, compute energy `H = 0.5*vpar^2 + mu*Bmod` and ensure matches reference value from prior run (hard-code expected value with tolerance 1e-10).
+     - Validate orthogonality: `dot_product(hcov, tracer%f%dhth)` should be near zero (tolerance 1e-12).
+   - Register in CMake and run: `ctest -R test_field_can_meiss_vmec`.
+
+3. **EQDSK guard (WILL_FAIL initially)** (`test/tests/test_field_can_meiss_eqdsk.f90`)
+   - Implement same checks as VMEC guard; currently expect failure.
+   - Mark test `WILL_FAIL TRUE` in CMake.
+   - Run: `ctest -R test_field_can_meiss_eqdsk` and confirm it fails but suite continues.
+
+4. **Geoflux field smoke** (`ctest -R test_magfie_geoflux`) – must pass.
+
+5. **Golden record awareness**
+   - Inspect `test/golden_record/canonical/simple.in` – note it uses `isw_field_type=0` (canonical flux).
+   - Plan to add new golden case for Meiss EQDSK later (see §4).
+
+---
+## 1. libneo Coordinate System (follow ../libneo/TODO.md)
+Complete libneo work first. Do not proceed until `../libneo/build/ctest -R test_coordinate_systems` passes and exposes `make_vmec_coordinate_system` / `make_geoflux_coordinate_system`.
+
+---
+## 2. Integrate coordinate_system_t (no wrappers)
+
+1. **Imports**
+   - Add `use libneo_coordinates` wherever geometry is needed: `src/field/field_can_meiss.f90`, `src/field_can.f90`, `app/simple_diag_meiss.f90`, `src/diag/diag_meiss.f90`, and any helper modules use the geometry.
+
+2. **Global storage** (`field_can_meiss`):
+   - Replace `field_noncan` with:
+     ```fortran
+     class(coordinate_system_t), allocatable :: meiss_coords
+     class(MagneticField), allocatable :: source_field
+     ```
+
+3. **Initialise coordinate system** (`init_meiss`)
+   - After storing incoming magnetic field, allocate coordinate system directly:
+     ```fortran
+     if (allocated(meiss_coords)) deallocate(meiss_coords)
+     if (geoflux_ready()) then
+       call make_geoflux_coordinate_system(meiss_coords)
+     else
+       call make_vmec_coordinate_system(meiss_coords)
+     end if
+     ```
+   - Note: `geoflux_ready()` should be a new public function from `field_geoflux`. Add it if missing (e.g., module procedure returning logical flag).
+
+4. **Use coordinate system** across the module:
+   - Replace all VMEC-specific data access (e.g., global arrays `lam_phi`, `h_cov`, etc.) with method calls:
+     - `meiss_coords%evaluate_point(u, position)`.
+     - `meiss_coords%covariant_basis(u, basis)` to compute derivatives.
+     - `meiss_coords%metric_tensor(u, g, ginv, sqrtg)`.
+   - Compute derivatives/invariants using the returned basis/metric—*no* direct `splint_vmec_data` or `geoflux_to_cyl` inside SIMPLE.
+
+5. **Magnetic field sampling**
+   - Continue using the existing `MagneticField` polymorphic object (`source_field`) for `evaluate` calls. Geometry now depends solely on `meiss_coords`.
+
+6. **Entry points update** (`src/field_can.f90`)
+   - Ensure `init_field_can` passes the magnetic field so `init_meiss` can allocate the right coordinate system.
+   - No new SIMPLE types are created; everything uses libneo’s `coordinate_system_t` instance.
+
+7. **Diagnostics**
+   - In `app/simple_diag_meiss.f90`, after reading config, call `field_from_file` and then `init_meiss(field)`; remove direct `init_vmec` + manual cache calls.
+   - `src/diag/diag_meiss.f90` simply relies on `rh_can`; confirm that `rh_can` internally uses `meiss_coords`.
+
+8. **Upgrade EQDSK test**
+   - Remove `WILL_FAIL TRUE` from EQDSK test once refactor is complete.
+   - Enhance checks in both VMEC and EQDSK tests to compare outputs against reference values stored in a small JSON/Fortran data file (e.g., `test/reference/meiss_expected.dat`). Populate this file by running the old VMEC pipeline and recording `Bmod`, `hth`, `hph` for a few canonical points; these become the expected values with tolerance.
+
+9. **Run focused unit tests**
+   ```bash
+   cmake --build build --target test_field_can_meiss_vmec.x test_field_can_meiss_eqdsk.x
+   ctest -R test_field_can_meiss_vmec
+   ctest -R test_field_can_meiss_eqdsk
+   ctest -R test_magfie_geoflux
+   ```
+   All must pass (<1s each).
+
+---
+## 3. Regression Tests & Golden Record
+1. **Add Meiss EQDSK golden test**
+   - Create `test/golden_record/meiss_eqdsk/simple.in` with EQDSK-specific configuration (e.g., `netcdffile='EQDSK_I.geqdsk'`, `isw_field_type=3` to request Meiss, set deterministic seeds).
+   - Update `test/golden_record/run_golden_tests.sh` so it copies the EQDSK file into each test case directory (download once to `$TEST_DATA_DIR/EQDSK_I.geqdsk`).
+   - Ensure symlink/hard copy placed in each run directory (similar to `wout.nc`).
+   - Add validation at the end of each run: check `diag_meiss_*` outputs exist if expected.
+
+2. **Golden record reference update**
+   - Generate reference outputs by running `test/golden_record/golden_record.sh main` (or chosen reference commit) after new code is in place. Keep `RUN_DIR_REF` results for diff.
+   - Ensure `compare_golden_results.sh` accounts for new files (update script to include Meiss-specific outputs).
+
+3. **System regression command**
+   - After every change, run:
+     ```bash
+     ./test/golden_record/golden_record.sh main
+     ```
+     and ensure the current run matches reference (no diffs aside from known improvements).
+
+---
+## 4. Diagnostics Automation
+1. **CTest driver** `test/diag/run_diag_meiss.cmake` (see prior template). Extend to check *both* VMEC and EQDSK outputs and to verify key numeric summaries (e.g., parse `diag_meiss_*.pdf` metadata using a small Python script or check associated data files if generated).
+2. **Register test** `diag_meiss_runs` with `TIMEOUT 120`.
+3. **Run** `ctest -R diag_meiss_runs`.
+
+---
+## 5. Examples & Documentation
+1. **Makefile** `examples/tokamak/Makefile` – ensure targets `all`, `run`, `diag`, `clean` as previously outlined.
+2. **README** `examples/tokamak/README.md` must include:
+   - Download step (`make`).
+   - Execution (`make run`), expected runtime, note about OpenMP threads.
+   - Diagnostics (`make diag`), location of generated plots.
+   - Troubleshooting tips (e.g., missing libneo branch).
+3. **Manual run**: follow README to completion; confirm outputs align with documentation.
+
+---
+## 6. Final Regression Sweep
+1. **Unit tests**: `ctest --output-on-failure` in SIMPLE repo.
+2. **System tests**: `./test/golden_record/golden_record.sh main` (ensure 0 exit code).
+3. **libneo tests**: `ninja -C ../libneo/build && ctest`.
+4. **Clean status**: `git status` (SIMPLE + libneo).
+5. **Document** results (changelog, PR summary).
+
+Completion criteria: unit & system tests all green, new golden record case passes, diagnostics produce EQDSK Meiss plots, documentation updated.

cmake/Util.cmake

@@ -1,5 +1,4 @@
 include(FetchContent)
-
 function(get_branch_or_main REPO_URL REMOTE_BRANCH)
     execute_process(
         COMMAND git rev-parse --abbrev-ref HEAD

examples/tokamak/.gitignore

@@ -0,0 +1 @@
+EQDSK_I.geqdsk

examples/tokamak/Makefile

@@ -0,0 +1,19 @@
+EQDSK := EQDSK_I.geqdsk
+URL   := https://crppwww.epfl.ch/~sauter/benchmark/EQDSK_I
+RUN   := ../../build/simple.x
+
+.PHONY: all run clean
+
+all: $(EQDSK)
+
+$(EQDSK):
+	@printf 'Downloading %s...\n' $(EQDSK)
+	@curl -L --fail --show-error --output $(EQDSK).tmp $(URL)
+	@mv $(EQDSK).tmp $(EQDSK)
+	@printf 'Saved %s\n' $(EQDSK)
+
+run: $(EQDSK)
+	@$(RUN) simple.in
+
+clean:
+	@rm -f $(EQDSK)

examples/tokamak/simple.in

@@ -0,0 +1,7 @@
+&config
+trace_time = 1d-2
+sbeg = 0.3d0
+ntestpart = 32
+netcdffile = 'EQDSK_I.geqdsk'
+isw_field_type = 1
+/

fpm.toml

@@ -6,7 +6,7 @@ maintainer = "albert@tugraz.at"
 copyright = "Copyright 2025, Plasma Physics at TU Graz"
 
 [dependencies]
-libneo.git = "https://github.com/itpplasma/libneo.git"
+libneo.git = { url = "https://github.com/itpplasma/libneo.git", branch = "tokamak" }
 openmp = "*"
 hdf5 = "*"
 netcdf = "*"

python/CMakeLists.txt

@@ -85,8 +85,10 @@ set(FILES_TO_WRAP
     magfie_wrapper.f90
 )
 
-# Add libneo files to wrap
-if(DEFINED ENV{CODE})
+# Add libneo files to wrap.
+# Prefer explicit CODE only when it is non-empty; otherwise use the resolved
+# libneo checkout from FetchContent (libneo_SOURCE_DIR).
+if(DEFINED ENV{CODE} AND NOT "$ENV{CODE}" STREQUAL "")
     set(LIBNEO_SOURCE_DIR $ENV{CODE}/libneo)
 else()
     set(LIBNEO_SOURCE_DIR ${libneo_SOURCE_DIR})

src/CMakeLists.txt

@@ -11,6 +11,7 @@ set(SOURCES
     field/field_coils.f90
     field/field_vmec.f90
     field/field_splined.f90
+    field/field_geoflux.f90
     field/vmec_field_eval.f90
     field/field_newton.F90
     field.F90
@@ -68,6 +69,7 @@ endif()
 
 add_library (simple STATIC ${SOURCES})
 
+
 # Link pyplot library to SIMPLE
 target_link_libraries(simple PUBLIC pyplot)
 
@@ -88,6 +90,7 @@ target_link_libraries(simple PUBLIC
 
 target_link_libraries(simple PUBLIC
     LIBNEO::neo
+    LIBNEO::magfie
 )
 
 # Conditionally link GVEC if enabled

src/coordinates/coordinates.f90

@@ -4,6 +4,9 @@ module simple_coordinates
     use vmec_coordinates, only: vmec_to_cyl_lib => vmec_to_cyl, &
                                  vmec_to_cart_lib => vmec_to_cart, &
                                  cyl_to_cart_lib  => cyl_to_cart
+    use geoflux_coordinates, only: geoflux_to_cyl_lib => geoflux_to_cyl, &
+                                     geoflux_to_cart_lib => geoflux_to_cart, &
+                                     cyl_to_geoflux_lib => cyl_to_geoflux
 
     implicit none
 
@@ -45,6 +48,33 @@ subroutine transform_cyl_to_cart(xfrom, xto, dxto_dxfrom)
 end subroutine transform_cyl_to_cart
 
 
+subroutine transform_geoflux_to_cyl(xfrom, xto, dxto_dxfrom)
+    real(dp), intent(in) :: xfrom(3)
+    real(dp), intent(out) :: xto(3)
+    real(dp), intent(out), optional :: dxto_dxfrom(3,3)
+
+    call geoflux_to_cyl_lib(xfrom, xto, dxto_dxfrom)
+end subroutine transform_geoflux_to_cyl
+
+
+subroutine transform_geoflux_to_cart(xfrom, xto, dxto_dxfrom)
+    real(dp), intent(in) :: xfrom(3)
+    real(dp), intent(out) :: xto(3)
+    real(dp), intent(out), optional :: dxto_dxfrom(3,3)
+
+    call geoflux_to_cart_lib(xfrom, xto, dxto_dxfrom)
+end subroutine transform_geoflux_to_cart
+
+
+subroutine transform_cyl_to_geoflux(xfrom, xto, dxto_dxfrom)
+    real(dp), intent(in) :: xfrom(3)
+    real(dp), intent(out) :: xto(3)
+    real(dp), intent(out), optional :: dxto_dxfrom(3,3)
+
+    call cyl_to_geoflux_lib(xfrom, xto, dxto_dxfrom)
+end subroutine transform_cyl_to_geoflux
+
+
 function get_transform(from, to)
     procedure(transform_i), pointer :: get_transform
     character(*), intent(in) :: from, to
@@ -56,6 +86,8 @@ function get_transform(from, to)
         select case (trim(to))
         case ('cart')
             get_transform => transform_cyl_to_cart
+        case ('geoflux')
+            get_transform => transform_cyl_to_geoflux
         case default
             call handle_transform_error(from, to)
         end select
@@ -68,6 +100,15 @@ function get_transform(from, to)
         case default
             call handle_transform_error(from, to)
         end select
+    case ('geoflux')
+        select case (trim(to))
+        case ('cyl')
+            get_transform => transform_geoflux_to_cyl
+        case ('cart')
+            get_transform => transform_geoflux_to_cart
+        case default
+            call handle_transform_error(from, to)
+        end select
     case default
         print *, "get_transform: Unknown transform from ", from
         error stop