Skip to content

Umbrella: replace avoidable numerical dependencies with AD-ready fortnum #291

@krystophny

Description

@krystophny

Coordination

This issue is the single entry point for the coordinator driving the
fortnum extraction. The coordinator dispatches each leaf issue per its
Model field, enforces the review ladder, updates tracker checklists, and
implements nothing directly. Coordinator model: fable (opus acceptable for
pure dispatch turns).

Fan-out

Tracker Repo Unblocked today? First actionable leaves
Phase 1: land the license/cleanup stack (#300) itpplasma/libneo yes #292 (sonnet), #293 (opus)
Phase 2/3: consume fortnum (#301) itpplasma/libneo no, after Phase 1 + fortnum M1 #296, #297, #298
M0 scaffold (lazy-fortran/fortnum#30) lazy-fortran/fortnum (private) yes leaves under #30
M1 math kit (lazy-fortran/fortnum#31) lazy-fortran/fortnum (private) no, after M0 leaves under #31
M2 ode (lazy-fortran/fortnum#32) lazy-fortran/fortnum (private) no, after M1 leaves under #32
M3 roots-interp-rng (lazy-fortran/fortnum#33) lazy-fortran/fortnum (private) no, after M1 leaves under #33
M4 integrate (lazy-fortran/fortnum#34) lazy-fortran/fortnum (private) no, after M1 leaves under #34
M5 hardening (lazy-fortran/fortnum#35) lazy-fortran/fortnum (private) no, last leaves under #35
Migrate avoidable numerics (NEO-2#88) itpplasma/NEO-2 no, after fortnum special/roots/integrate spawn on milestone close
Migrate avoidable numerics (KAMEL#134) itpplasma/KAMEL no, after fortnum special/quadrature/ode spawn on milestone close
Drop GSL/FFTW copies (MEPHIT#24) itpplasma/MEPHIT no, after fortnum interp/ode spawn on milestone close
Migrate generic numerics (GORILLA#43) itpplasma/GORILLA no, after fortnum ode spawn on milestone close
Adopt fortnum RNG/ODE helpers (SIMPLE#371) itpplasma/SIMPLE no, after fortnum rng/ode spawn on milestone close
Replace vendored helper copies (NEO-RT#43) itpplasma/NEO-RT no, after fortnum interp/ode spawn on milestone close

fortnum is private during incubation. Reference its issues by bare number
only; quote no private content into public-repo bodies.

Execution order

  1. libneo Phase 1 (Tracker: Phase 1 — land the license/cleanup stack #300) and fortnum M0 (Test even order in spline three to five #30) run in parallel today.
  2. fortnum M1 after M0.
  3. libneo Phase 2 (Tracker: Phase 2/3 — consume fortnum #301, leaves Consume AD-ready fortnum via CMake FetchContent #296 and Replace src/math with AD-ready fortnum + keep neo_* compatibility wrappers #297) after Phase 1 and M1.
  4. fortnum M2-M4 and libneo Wrap odeint_allroutines around AD-ready fortnum_ode #298.
  5. Downstream trackers spawn leaves as their fortnum milestones close.
  6. fortnum M5 and libneo Remove compatibility wrappers after downstream migration #299 last.

Execution conventions

Model ladder (implementation):

  • haiku: label/milestone housekeeping, checklist updates, re-export boilerplate. No algorithms.
  • sonnet: fully specified ports (source path + SHA given), wrappers, CI fixes with a clear failing log.
  • opus: build-system and CI scaffolding, multi-file refactors, ports that change APIs.
  • fable: API design issues, clean-room algorithm work, review of opus-implemented PRs.

Review ladder: sonnet PRs reviewed by opus; opus PRs reviewed by fable; haiku changes gated by CI only.

PR rules:

  • One leaf issue = one PR. Branch from the base named in the issue.
  • PR title = issue title. Body: Closes #N plus a ## Verification section with real failing-before and passing-after output.
  • Squash merge. Stacked PRs only where the issue names a non-main base.

Executor rule: if any issue under "Blocked by" is open, stop and report instead of starting.

Standing instructions for the coordinator

  • If a leaf's Blocked-by is open, skip it.
  • Create downstream leaf issues from the tracker tables only when the named
    fortnum milestone closes. Use the leaf template below.
  • Never link private fortnum content into public-repo issue bodies beyond
    bare issue references.
  • Keep the Fan-out table's "Unblocked today?" column current as milestones
    close.
Leaf issue template
Model: <haiku|sonnet|opus|fable>
Review: <per ladder>
Blocked by: <#N, ... or "none">
Branch from: <base>
Files: <exact paths>

## Spec
...

## Port from / Context
...

## Do not
...

## Verify
<exact commands and expected outcomes>

## Done when
<checklist>

Goal

Create lazy-fortran/fortnum: a small MIT-licensed numerical-computing library for generic algorithms that are currently duplicated, vendored, or supplied through avoidable upstream numerical dependencies in libneo and downstream plasma codes.

This is not a plan to replace the whole scientific software stack. We explicitly keep the heavy infrastructure dependencies upstream:

keep upstream:
  NetCDF / NetCDF-Fortran
  HDF5
  BLAS
  LAPACK
  MPI
  OpenMP
  Fortran compiler/runtime
  CMake
  Python build/packaging tools where needed

The target is narrower and realistic:

replace / absorb into fortnum:
  GSL / FGSL wrappers
  direct gsl_* C bindings
  FFTW usage where local clean-room FFT is good enough
  old LGPL/Burkardt/IQPACK-style quadrature generators
  old QUADPACK/SLATEC/AMOS/libcerf vendoring where it is uncontrolled or unnecessary
  duplicated RKF45 / Cash-Karp / odeint code
  root-finding wrappers
  generic integration/quadrature helpers
  generic interpolation/polynomial/search helpers
  generic RNG/sampling wrappers

libneo should remain the plasma/domain library: VMEC/Boozer/GEQDSK, magnetic coordinates, magfie, species/collisions/transport domain interfaces, Python domain I/O, and compatibility wrappers. fortnum should contain generic numerical kernels reusable by SIMPLE, NEO-2, KAMEL, GORILLA, and lazy-fortran packages.

Proposed tagline:

fortnum: MIT-licensed numerical computing for modern Fortran.

Why

Recent cleanup work already points in this direction:

The same pattern exists across downstream codes: generic ODE, quadrature, special-function, root-finding, interpolation, and RNG code is duplicated, vendored, or tied to GSL/FFTW/FGSL interfaces.

The point is supply-chain control for the avoidable numerical layer, not replacing foundational ecosystem libraries like HDF5, NetCDF, BLAS, or LAPACK.

External design references

These are the design points worth copying, not necessarily their exact APIs:

  • SciPy API structure: public submodules are explicitly defined, stable, and deprecations are used before incompatible changes. Submodules should be as self-contained as possible with minimal inter-submodule dependencies.
    https://docs.scipy.org/doc/scipy/reference/index.html
  • SciPy thread-safety guidance: thread workers should preferably own their own arrays/objects; functions should not mutate user arrays unless explicitly documented; stateful objects require clearer concurrency semantics.
    https://docs.scipy.org/doc/scipy/tutorial/thread_safety.html
  • SUNDIALS design direction: solver objects, vector abstractions, pluggable linear/nonlinear solvers, modern Fortran interfaces, and hidden implementation changes without breaking user APIs.
    https://arxiv.org/abs/2011.10073
  • SUNDIALS performance direction: fused vector operations, local/global reduction separation, and backend-flexible data structures for parallel/heterogeneous systems.
    https://arxiv.org/abs/1909.12966
  • GSL scope reference: special functions, integration, FFT, roots, minimization, ODEs, interpolation, RNG, statistics, splines, sparse matrices, etc. GSL itself is GPL and therefore should be an optional oracle only, not a shipped dependency.
    https://www.gnu.org/software/gsl/
  • FFTW license reference: GPL by default, with non-free alternative licensing; therefore also oracle-only or replaced by clean-room/permissive code in shipped artifacts.
    https://www.fftw.org/doc/License-and-Copyright.html
  • QUADPACK reference: public-domain Fortran quadrature family and useful semantic/API model for adaptive integration.
    https://www.netlib.org/quadpack/
  • PocketFFT reference: permissively licensed FFT implementation used in the Python ecosystem; useful as an algorithm/design comparison, but do not copy unless license/provenance is explicitly checked.
    https://github.com/mreineck/pocketfft

Naming decision

Use fortnum, not fortmath, fortcalc, fortalg, or fortsl.

Rationale:

  • fortnum means numerical computing, which matches the scope.
  • fortmath is broader and vaguer: symbolic math, algebra, elementary math, discrete math.
  • fortcalc sounds like a calculator/calculus package and is too narrow.
  • fortalg is ambiguous between algorithms and algebra.
  • fortsl is too close to GSL and not descriptive enough.

Public module convention:

use fortnum_special
use fortnum_integrate
use fortnum_quadrature
use fortnum_fft
use fortnum_ode
use fortnum_roots
use fortnum_interp
use fortnum_rng
use fortnum_linalg_lite

Private implementation modules should use a clear private prefix, for example:

use fortnum_ode_cash_karp_impl
use fortnum_fft_bluestein_impl

Only documented top-level modules are public API.

What stays upstream and is not part of this issue

Do not replace these in this effort:

NetCDF / NetCDF-Fortran
HDF5
BLAS
LAPACK
MPI
OpenMP
Fortran compiler/runtime
CMake
Python packaging/build tooling

Do not move these to fortnum:

VMEC/Boozer/GEQDSK/Hamada coordinate logic
magfie physics interfaces
HDF5/NetCDF/Python domain I/O
plasma species/collision/transport domain models
orbit/tetrahedron/symplectic physics algorithms, except clearly generic kernels

Optional oracle/test-only dependencies are allowed:

GSL
FGSL
FFTW
SciPy
mpmath

but they must not be linked into shipped MIT artifacts.

Inventory: replace or absorb from libneo

A. Clean-room math kit from #290

Move after #286/#289/#290 settle, preserving oracle tests and benchmarks:

src/math/bessel_i.f90         -> fortnum_special
src/math/bessel_k.f90         -> fortnum_special
src/math/dawson.f90           -> fortnum_special
src/math/gauss_kronrod.f90    -> fortnum_integrate / fortnum_quadrature
src/math/gauss_quadrature.f90 -> fortnum_quadrature
src/math/fft.f90              -> fortnum_fft
lower incomplete gamma (#286) -> fortnum_special_gamma

Keep temporary libneo wrappers so existing code can still import old module names while downstream migration happens.

B. ODE integrators

Move generic ODE machinery:

src/odeint/odeint_allroutines.f90
test/odeint/*
doc/UserDoc/ODE_Integration.tex material that is generic rather than libneo-specific

Current code is generic adaptive Cash-Karp RK5(4) and event detection. This belongs in fortnum_ode, but not as a raw copy. Refactor to a clean API first:

type(ode_cash_karp_t) :: solver
type(ode_status_t) :: status

call solver%init(nvar)
call solver%integrate(y, x1, x2, rtol, atol, derivs, status)

Also provide a procedural wrapper:

call odeint_cash_karp(y, x1, x2, rtol, derivs, status)

Keep in libneo:

module odeint_allroutines_sub
  use fortnum_ode
  ! old argument order and old procedure names
end module

C. RNG utilities

Candidate:

src/MC/rng.f90

Do not merely move the existing wrapper around random_number. Design fortnum_rng around explicit RNG state so parallel runs are reproducible and thread-safe:

type(rng_t) :: rng
call rng%seed(seed, stream_id)
call rng%uniform(u)
call rng%normal(x)

Initial compatibility wrapper can preserve getran(irand, ur).

D. Interpolation / polynomial / search helpers

Candidate generic pieces:

src/plag_coeff.f90
src/binsrc.f90
selected generic interpolation kernels under src/interpolate/ if not VMEC/magfie/domain-specific

Do not move domain-specific spline construction that knows about VMEC/Boozer/magnetic coordinates. Move only generic interpolation, search, and Lagrange/polynomial kernels.

E. Root finding / optimization placeholders

libneo itself currently does not seem to own the main root-finding wrappers, but downstream codes do. fortnum_roots should be created early enough to absorb NEO-2 and KAMEL GSL wrappers.

Minimum public API:

call root_bisection(f, a, b, xtol, rtol, root, status)
call root_newton(fdf, x0, xtol, rtol, root, status)
call root_brent(f, a, b, xtol, rtol, root, status)  ! later

No global procedure pointers; pass callbacks through explicit interfaces and optional context objects.

Inventory: replace or absorb from downstream codes

SIMPLE

Search hits indicate current/local use of:

src/samplers.f90
src/collis_alphas.f90 random-number use
orbit/canonical-coordinate RK-related code paths and tests

Action: after fortnum_rng and fortnum_ode exist, audit SIMPLE samplers and any generic ODE/RK helpers for migration. Keep symplectic/orbit integrators in SIMPLE if they are physics-specific.

NEO-2

Replace avoidable GSL/FGSL dependencies:

COMMON/gsl_specialfunctions_mod.f90   -> fortnum_special
COMMON/gsl_roots_routines_mod.f90     -> fortnum_roots
COMMON/gsl_integration_routines_mod.f90 -> fortnum_integrate
COMMON/gsl_bspline_routines_mod.f90   -> fortnum_interp / fortnum_bspline
cmake/fgsl.cmake                      -> remove once no longer needed

Also audit:

COMMON/rk4_kin.f90
COMMON/rk4_kin_mod.f90
COMMON/flint_mod.f90
NEO-2-QL/flint.f90
NEO-2-PAR/flint.f90

for generic ODE/integration pieces versus domain-specific kinetic-field logic.

Target: replace FGSL/GSL dependency where possible with fortnum_special, fortnum_integrate, fortnum_roots, and fortnum_interp.

KAMEL

Replace avoidable GSL/direct C bindings:

KIM/src/math/gsl_mod.f90     -> fortnum_special
cmake/FetchGSL.cmake         -> remove once no longer needed

Audit/replace quadrature and ODE code:

KIM/src/math/quadpack_integration_m.f90
KIM/src/math/quadpack_compat_xerror.f90
KIM/src/kernels/integrals.f90
KIM/src/kernels/integrals_adaptive.f90
KIM/src/kernels/integrands_adaptive.f90
KIM/src/math/ddeabm/*
KiLCA/math/odeint_allroutines.f
KiLCA/solver/VER_5_STABLE/odeint_allroutines.f
KiLCA/flre/conductivity/odeint_allroutines.f
common/equil/rk4_integrator.f90
PreProc/fourier/src/rk4dc.f90

Audit special-function vendoring/provenance:

KiLCA/math/bessel/amos/*
KiLCA/math/bessel/644/*
common/math/slatec/*
KIM/src/math/libcerf-main/*
KIM/src/math/use_libcerf_mod.f90
KiLCA/math/hyper/hyper1F1.cpp

These should become either:

A. replaced by clean-room fortnum_special routines, or
B. deliberately vendored with explicit license/provenance and tests, if acceptable.

Prefer A for the long-term supply-chain-control goal.

Audit duplicate interpolation helpers:

KIM/src/util/plag_coeff.f90
QL-Balance/src/base/plag_coeff.f90
KiLCA/solver/VER_5_STABLE/plag_coeff.f90

Target: domain integrands stay in KAMEL; generic quadrature/ODE/special-function kernels move to fortnum.

GORILLA

Candidates:

SRC/odeint_rkf45.f90
SRC/runge_kutta_mod.f90
SRC/pusher_tetra_rk.f90
SRC/plag_coeff.f90

Target: use fortnum_ode for RKF45/Cash-Karp-like utilities where practical; keep orbit/tetrahedron physics in GORILLA.

MEPHIT / NEO-RT / benchmark_orbit / sparse_draft

From #258:

MEPHIT: trivial FFTW interface copy and diverged/orphaned field-line integration copies
NEO-RT: historical copies/symlinks of binsrc, plag_coeff, spline helpers
benchmark_orbit: pre-built archives
sparse_draft: broken raw download of hdf5_tools

Action: after fortnum_interp and fortnum_ode exist, audit these for use of binsrc, plag_coeff, spline helpers, and old ODE code. Do not move hdf5/domain I/O into fortnum.

fortnum package structure

Initial repository layout:

fortnum/
  CMakeLists.txt
  fpm.toml                    # optional only if workable; CMake remains primary for subtargets
  src/
    fortnum_kinds.f90
    fortnum_status.f90
    special/
      fortnum_special.f90
      fortnum_special_gamma.f90
      fortnum_special_bessel.f90
      fortnum_special_error.f90
      fortnum_special_dawson.f90
      fortnum_special_faddeeva.f90
    quadrature/
      fortnum_quadrature.f90
      fortnum_integrate.f90
      fortnum_integrate_gk.f90
    fft/
      fortnum_fft.f90
    ode/
      fortnum_ode.f90
      fortnum_ode_cash_karp.f90
      fortnum_ode_events.f90
    roots/
      fortnum_roots.f90
    interp/
      fortnum_interp.f90
      fortnum_polynomial.f90
      fortnum_bspline.f90
    rng/
      fortnum_rng.f90
    linalg_lite/
      fortnum_linalg_lite.f90
  test/
    special/
    quadrature/
    fft/
    ode/
    roots/
    interp/
    rng/
    oracle/                  # optional GSL/FFTW/SciPy/mpmath comparisons
  benchmark/
  docs/
    api.md
    design.md
    migration_libneo.md

Interface rules

Public API

  • Public modules are exactly the documented top-level fortnum_* modules.
  • Implementation modules must be private by convention and not imported by downstream codes.
  • Deprecations before incompatible changes once downstream codes use it.
  • Avoid exposing algorithm internals unless they are stable concepts, e.g. ode_cash_karp_t is public, but rk_stage_buffer_t should not be.

Thread safety

  • No hidden global mutable state in public algorithms.
  • No module-level save state for solver status, callbacks, work arrays, RNG seeds, FFT plans, or error flags.
  • Use caller-owned workspaces/solver objects for stateful algorithms.
  • Pure/elemental scalar special functions where possible.
  • Caller-owned, read-only plan objects for FFT.
  • Explicit RNG state per thread/stream.
  • Parallel tests must cover simultaneous calls from OpenMP regions.
  • If an object is not safe for concurrent mutation, document it and detect/reject illegal concurrent calls where feasible.

Performance

  • No allocation in scalar hot paths.
  • Work arrays are preallocated in solver/workspace objects.
  • Provide array/batched APIs for special functions and FFT when it avoids repeated setup overhead.
  • Keep inner loops simple and vectorizable: contiguous arrays, do concurrent/!$omp simd where justified by tests.
  • Avoid polymorphism in inner loops; use generic interfaces at API boundaries, then monomorphic kernels internally.
  • Separate scalar, vector, and batched APIs where that improves performance clarity.
  • Benchmark against current libneo/GSL/FFTW behavior before and after migration.

Error handling

  • No error stop in reusable numerical kernels except for impossible internal invariants.
  • Return explicit status objects or integer status codes:
type :: fortnum_status_t
  integer :: code
  character(len=:), allocatable :: message
  integer :: n_eval
  integer :: n_iter
  real(dp) :: estimated_error
end type
  • Fatal behavior can be a wrapper choice, not the kernel default.

Licensing/provenance

  • src/ must be MIT-compatible only.
  • GSL/FFTW/FGSL are optional oracle-test dependencies only.
  • Do not copy Numerical Recipes code.
  • For AMOS/SLATEC/QUADPACK/PocketFFT/libcerf-derived material, audit provenance and license explicitly before moving.
  • Prefer clean-room implementation for long-term control.
  • Generated coefficient tables must include a reproducible script and reference source.

Migration plan

Phase 1: stabilize libneo-side cleanup

Phase 2: create fortnum and move clean-room math

  • Create lazy-fortran/fortnum under MIT.
  • Move clean-room special functions/quadrature/FFT from Add clean-room math kit; drop FFTW link and LGPL quadrature #290.
  • Add old neo_* compatibility wrappers in libneo.
  • Add oracle tests and benchmarks in fortnum.
  • Make libneo consume fortnum through CMake FetchContent or another controlled dependency mechanism.

Phase 3: move ODE machinery

  • Design fortnum_ode object/workspace API.
  • Port Cash-Karp RK5(4), adaptive controller, event detection, and tests.
  • Keep libneo odeint_allroutines wrapper.
  • Replace GORILLA/SIMPLE/KAMEL generic RKF45/ODE copies where feasible.

Phase 4: roots/integration/interpolation/RNG

  • Add fortnum_roots bisection/Newton/Brent-style APIs.
  • Replace NEO-2 FGSL root wrappers.
  • Extend fortnum_integrate toward QAG/QAGS/QAGP/QAGIU-compatible semantics.
  • Replace NEO-2 FGSL integration wrappers where feasible.
  • Add fortnum_interp/fortnum_polynomial; migrate plag_coeff, binsrc, and duplicate polynomial helpers.
  • Add fortnum_rng; migrate trivial RNG wrappers and establish reproducibility policy.

Phase 5: special functions beyond current libneo needs

  • Replace KAMEL direct GSL bindings for bessel_In, erf, erfc, dawson.
  • Audit and replace or explicitly vendor KAMEL SLATEC/AMOS/libcerf special-function code.
  • Decide whether fortnum_special should own Faddeeva/plasma-dispersion-function functionality.

Phase 6: downstream gates and deprecation

  • Use Per-PR reverse-dependency gate against downstream codes #278-style per-PR downstream gate for libneo+fortnum changes.
  • Add fast downstream CI modes for SIMPLE, NEO-2, KAMEL, GORILLA.
  • Add nightly slow/golden/performance downstream jobs.
  • Deprecate libneo generic numerical module names only after downstream replacements exist.
  • Remove compatibility wrappers after at least one release cycle and green downstream CI.

Acceptance criteria

  • libneo can build without installed GSL and FFTW.
  • fortnum has no GPL/LGPL code in shipped source or linked artifacts.
  • GSL/FGSL/FFTW/SciPy/mpmath are optional oracle-test dependencies only.
  • NetCDF, HDF5, BLAS/LAPACK, MPI/OpenMP remain accepted upstream dependencies and are not in scope for replacement.
  • fortnum APIs are documented, stable, and separated into public/private modules.
  • ODE, FFT, quadrature, special-function, roots, interpolation, and RNG APIs are thread-safe by construction or explicitly documented otherwise.
  • No hidden global mutable state in reusable numerical kernels.
  • Benchmarks show acceptable regression bounds compared with current libneo/GSL/FFTW-backed paths.
  • SIMPLE, NEO-2, KAMEL, and GORILLA fast CI pass against libneo using fortnum.
  • Existing numerical regression/golden tests remain green or have documented, reviewed tolerance changes.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions