Skip to content

[spectre-02] Evaluate SPECTRE vector potential with analytic second derivatives #369

Description

@krystophny

Stage: libneo / field evaluation
Source language: Fortran 2008
Manuals to read first: SPECTRE fortran_src/base_functions_mod.F90 (get_cheby, get_zernike, get_zernike_d2 — recurrences, normalization by 1/(n+1), regularization); SPECTRE fortran_src/magnetic_field_mod.F90:32-115 (compute_magnetic_field — how TT scaling with sbar=(1+s)/2 and the half derivative factor enter); SPECTRE python spectre/postprocess/_get_field.py (reference implementation with derivatives); itpplasma/SPOIC spoic.f90 lines 35-100 (working Fortran evaluation loop for SPEC).
Depends on: #368

Goal

Given spectre_data_t, evaluate at a point (lvol, s, theta, zeta) the covariant vector potential components A_theta, A_zeta and all first and second derivatives with respect to (s, theta, zeta), analytically, for both interior volumes (Chebyshev radial basis) and the axis volume (Zernike radial basis with SPEC regularization factors). This is the only field-data evaluation primitive; everything downstream (B, metric coupling, guiding-center derivatives) consumes it.

Files to edit

  • src/spectre/spectre_basis.f90: NEW. Module spectre_basis: get_cheby_d2(s, lrad, T), get_zernike_d2(sbar, lrad, mpol, Z), eval_spectre_vector_potential(data, lvol, s, theta, zeta, av) filling a result type with Ath, Azt, dAth(3), dAzt(3), d2Ath(6), d2Azt(6) (order: ss, st, sz, tt, tz, zz).
  • test/spectre/test_spectre_basis.f90: NEW.
  • test/spectre/data/reference_vecpot.dat: NEW. Reference values exported from SPECTRE's python get_vec_pot/_get_field.py for the committed fixture at ~20 sample points per volume; generation added to tools/generate_spectre_test_file.py.

Behavior to implement

  • Fourier sum over mn modes with arg = im(ii)*theta - in(ii)*zeta; even (Ate/Aze with cos) and odd (Ato/Azo with sin) parts.
  • Interior volumes: Chebyshev recurrence in s on [-1,1] with value/first/second derivative, SPEC normalization conventions applied exactly as in base_functions_mod.F90.
  • Axis volume (lvol == 1): Zernike in sbar = (1+s)/2 indexed by radial degree and poloidal m; chain-rule factors from d sbar/ds = 1/2; the m-dependent regularization keeps all outputs finite and correct as s -> -1.
  • Everything pure/threadsafe: no module state; OpenMP-hot-loop safe.

Scaffold

do ii = 1, data%mn
    arg = data%im(ii)*theta - data%in(ii)*zeta
    carg = cos(arg); sarg = sin(arg)
    ! radial factor T(0:2) from cheby or zernike(:, im(ii), 0:2)
    do ll = 0, data%Lrad(lvol)
        av%Ath = av%Ath + (Ate(ll,ii)*carg + Ato(ll,ii)*sarg)*T(ll,0)
        av%dAth(1) = av%dAth(1) + (...)*T(ll,1)   ! etc. for all 8 derivative slots
    end do
end do

Acceptance scenarios (BDD)

  • Given the fixture and the reference sample points, when evaluating A_theta, A_zeta, then values match the SPECTRE python reference to 1e-12 in every volume including the axis volume.
  • Given central finite differences of the evaluated A over a step sweep h in {1e-3..1e-6}, then analytic first and second derivatives converge with the expected order (verifies every one of the 18 derivative slots).
  • Given points approaching the axis (s -> -1 in volume 1), then all outputs remain finite and FD-consistent (Zernike regularization correct).
  • Given |s| slightly > 1 (Newton overshoot), then evaluation returns the smooth polynomial extension without error (documented behavior, needed by event location downstream).

Success criteria

make test 2>&1 | grep -E 'spectre_basis.*(Passed|OK)'

Non-goals

  • No metric/geometry (spectre-03), no B or Bmod (spectre-04).
  • No performance work (no batching, no splines) — correctness first; hot-loop layout may be revisited later.

Verification

cmake --build build && make test

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Fields

    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