batch_interpolate_3d: kink-robust BC variant + compact coefficient storage

[krystophny]

batch_interpolate_3d: kink-robust BC variant + compact coefficient storage (for real-space field splines)

Two enhancements to the batch 3D splines (src/interpolate/batch_interpolate_3d.f90)
motivated by adding an ASCOT5-style direct real-space (R,phi,z) magnetic-field spline as
a fast, GPU-friendly field representation for full-orbit tracers (itpplasma/SIMPLE#423;
perf context SIMPLE#422). The existing engine already covers most needs -- cubic AND
quintic orders, periodic/natural BCs, value/1st/2nd-derivative evaluate, SoA batch, and
GPU !$acc routine seq + resident device kernel with the managed-memory fix. Two gaps:

1. Kink-robust / boundary-aware spline variant

Today the only BCs are natural and periodic (spl_reg/spl_per in
spl_three_to_five.f90); there is no monotone/Steffen-limited/WENO/one-sided option. A
real-space field grid is rectangular but the plasma boundary (LCFS) is not, so a grid
line crossing rho=1 sees a C1 kink (real field inside, clamped/extrapolated outside). A
global C2 cubic rings near that kink within ~3 nodes; a quintic rings worse and wider --
and the kink sits exactly at the edge where particle losses are decided. Request: a
limited/one-sided (e.g. Steffen-monotone or boundary ghost-node) construction option,
selectable per axis, that suppresses the overshoot near a data discontinuity while
keeping the interior accuracy. Lets us spline up to the LCFS without corrupting the edge
shell, and is reusable for any field with a sharp boundary.

2. Compact coefficient storage option

The current batch layout stores a dense (order+1)^3 coefficient block per cell
(coeff(iq, k1,k2,k3, i1,i2,i3)). For a real-space field grid this is the dominant
memory cost: a 4-quantity QUINTIC field at 200x64x200 is ~190 GB; cubic ~36 GB. ASCOT5's
interp3Dcomp stores one compact set per node (function value + 2nd derivatives, 8 per
node for tricubic) and reconstructs the cell polynomial on the fly, which is far lighter
and is what makes high-resolution 3D field grids affordable. Request: an optional compact
storage mode for BatchSplineData3D (store f + second derivatives, reconstruct in
evaluate), so high-resolution and/or quintic real-space field splines fit in memory. The
device evaluate kernels already exist; the change is the storage/reconstruction path.

Both are opt-in and do not change the default natural/periodic dense path. Needed for the
SIMPLE real-space-field-spline work (SIMPLE#423) but generally useful for any large 3D
field spline. Related: the chartmap covariant-basis work in PR #325.

[krystophny]

Automated triage pass (Claude), grounded in the repository at fe3f76e. Verify before acting.

Still open and still accurate: neither feature exists at fe3f76e. But the issue bundles two independent work items with very different risk profiles, and part 1 as literally worded ("boundary ghost-node", "one-sided BC") does not address the problem it describes. Recommendation below: split into two atomic issues, and reframe part 1.

Purpose

Both requests come from the same driver, itpplasma/SIMPLE#423: a direct real-space (R, phi, z) magnetic-field spline for the Cartesian full-orbit pusher, so the hot loop stops paying a damped-Newton chartmap inversion per step (itpplasma/SIMPLE#422). That representation is only viable if (a) a high-resolution 3D grid fits in memory — on a GPU, where it must be device-resident — and (b) the spline does not ring where the rectangular box crosses the LCFS, because that ring sits exactly in the shell where particle losses are decided. ASCOT5 makes the same representation affordable with compact per-node storage; libneo has no equivalent.

Current state

Storage. BatchSplineData3D (src/interpolate/batch_interpolate_types.f90:43-58) has exactly one coefficient array:

      ! Memory layout: (n_quantities, 0:order1, 0:order2, 0:order3, n1, n2, n3) for
      ! cache efficiency
      real(dp), dimension(:, :, :, :, :, :, :), allocatable :: coeff

allocated dense at four sites in src/interpolate/batch_interpolate_3d.f90 — lines 131-132, 311-312, 697-698 and 718-719 — one per builder:

            allocate (spl%coeff(n_quantities, 0:N1_order, 0:N2_order, 0:N3_order, &
                                n1, n2, n3), stat=istat)

Evaluation is a 3D Horner sweep straight off that array (src/interpolate/spline3d_many_point_body.inc:41-74), included by evaluate_batch_splines_3d_many at batch_interpolate_3d.f90:971-973 and by the resident/_der/_der2 variants.

The issue's memory figures are wrong; the real ones still make the case. For nq = 4 on 200 x 64 x 200 = 2 560 000 nodes at 8 B:

	doubles per node per quantity	total
dense cubic, (3+1)^3	64	5.24 GB
dense quintic, (5+1)^3	216	17.69 GB
compact cubic, 2^3	8	0.66 GB
compact quintic, 3^3	27	2.21 GB

i.e. a clean 8x reduction in both cases, not the ~190 GB / ~36 GB quoted. 17.7 GB still does not fit an 80 GB device alongside the rest of a SIMPLE run at higher resolution, and the peak-memory work in 04e0e22 ("interpolate: build batch spline 3D coeffs in place to halve construction peak memory (#341)") shows this is already a live constraint.

Boundary conditions. The only options are natural and periodic. Confirmed by reading spl_three_reg_line (src/interpolate/acc_spline_build_order5.f90:60-96): it sets b(1)=d(1)=0 and c(n)=0, i.e. zero second derivative at both ends — natural, not not-a-knot. rg -in 'steffen|monotone|limiter|one_sided|fritsch|akima' src/ returns nothing.

There are three builders, each with its own BC dispatch, and all three must be touched by either feature:

builder	lines	BC dispatch
construct_batch_splines_3d_legacy	64-205	spl_reg/spl_per (src/spl_three_to_five.f90:587,639) at lines 152/154, 172/174, 194/196
construct_batch_splines_3d_lines	207-530	inline spl_{three,four,five}_{reg,per}_line calls at lines 347-497
construct_batch_splines_3d_resident_device_impl	582-885	spl_build_line_inplace (acc_spline_build_order5.f90:17-58, dispatches solely on order and periodic) at lines 802, 827, 852

construct_batch_splines_3d (line 45) is a thin selector: #ifdef _OPENACC it calls the resident-device path, otherwise _legacy. So the CPU default and the GPU default are different builders. _lines is separately public (line 22).

The reframing. periodic(3) and the natural end conditions act at the ends of each grid line. The LCFS crossing in a rectangular (R, phi, z) box is an interior data kink, typically many nodes from either end. A ghost-node or one-sided end condition therefore cannot suppress it. Only a local, data-dependent limiter (Steffen/Fritsch–Carlson-type) does, and the trade is larger than the issue admits:

• a Steffen-limited cubic Hermite is C1, not C2 — evaluate_batch_splines_3d_der2 along such an axis becomes discontinuous at nodes, which matters directly for SIMPLE#423 option (B), where grad|B| is read from _der2;
• it is third-order, not fourth-order, even where the data is smooth and the limiter never clips. Measured below: the sup-norm error ratio under halving is 8.0, not 16.

So the issue's "while keeping the interior accuracy" is not achievable. Issue B must state the loss of one order and of C2 explicitly.

Reproduction / evidence

Memory. No measurement is needed — size(spl%coeff) is exactly nq*(o1+1)*(o2+1)*(o3+1)*n1*n2*n3, so the numbers above are arithmetic, and the check belongs in a test as an assertion rather than in a scratch program. For the benchmark case used in the acceptance criteria, 120 x 32 x 120 (460 800 nodes), nq = 4, order = [5,5,5]:

• dense: 4 * 216 * 460800 = 398 131 200 doubles = 3.19 GB
• compact: 4 * 27 * 460800 = 49 766 400 doubles = 0.40 GB

Build and run the existing suite with:

cmake --preset default && cmake --build --preset default && cd build && ctest --output-on-failure

Kink. No fixture exists today. The 1D slice of the failure mode, verified numerically before writing this:

# 1D slice of the failure mode, to be ported to a Fortran test
import numpy as np
x = np.linspace(0.0, 1.0, 65)
xk = 0.35
f = np.where(x < xk, np.cos(2*np.pi*x), np.cos(2*np.pi*xk))   # slope jump at x = xk

The clamp point is 0.35, deliberately not an extremum of cos(2*pi*x). At x = 0.5 the derivative of cos(2*pi*x) already vanishes, so clamping there produces only a second-derivative jump and the natural cubic undershoots by a negligible 2.4e-4; that would be a near-useless regression guard. At x = 0.35 the slope genuinely jumps.

Measured (scipy CubicSpline(bc_type='natural'), evaluated on 8001 uniform points):

	data min	spline min	undershoot
natural cubic	-0.587785	-0.589363	1.58e-3, i.e. 0.10% of the data range
Steffen-limited Hermite	-0.587785	-0.587785	0.0 to round-off

So the ring is real and the limiter removes it completely, but the amplitude is ~0.1% of the range for a cubic, not the "2-4%" that an earlier draft of this comment asserted. Quintic is expected to ring wider and harder; that has not been measured here and should be measured in the PR rather than assumed.

Order of accuracy under the limiter, measured on the same reference Steffen implementation, sup norm over 1000 interior points:

n	f = exp(x) error	ratio	f = tanh(2x) error	ratio
33	1.083e-06		3.688e-06
65	1.364e-07	7.94	4.911e-07	7.51
129	1.702e-08	8.01	5.954e-08	8.25
257	2.127e-09	8.00	7.547e-09	7.89

Ratio 8 = O(h^3). Both test functions are strictly monotone on [0,1], so the limiter never clips; this is the best case. On a non-monotone smooth function (sin(3x)*exp(-x), which has an interior maximum) the same code gave ratios of 3-5, because the limiter clips at the extremum. The convergence test must therefore use a monotone function, or it measures the limiter's clipping, not its accuracy.

Behaviour (BDD)

Scenario: compact storage reproduces the dense spline to round-off
  Given a random 4-quantity dataset on a 33 x 17 x 33 grid
  And spl_dense built by construct_batch_splines_3d with order=[3,3,3], periodic=[F,T,F]
  And spl_compact built with the same arguments plus compact=.true.
  When both are evaluated at 10000 pseudo-random interior points
  Then max|y_compact - y_dense| <= 1e-12 * max|y_dense|
  And size(spl_compact%coeff_compact) == 8 * 4 * 33 * 17 * 33
  And spl_compact%coeff is not allocated

Scenario: compact storage supports quintic order
  Given the same dataset built with order=[5,5,5] and compact=.true.
  When evaluated with evaluate_batch_splines_3d_many_der at 10000 points
  Then value and all three first derivatives match the dense build to 1e-11 relative
  And size(spl_compact%coeff_compact) == 27 * 4 * 33 * 17 * 33

Scenario: Steffen-limited axis does not overshoot a clamped plateau
  Given y(i) = cos(2*pi*x(i)) for x < 0.35 and y = cos(2*pi*0.35) for x >= 0.35
  And 65 uniform nodes on [0,1]
  And order(1)=3 with bc(1)=BC_STEFFEN, other axes natural
  When evaluated on 8001 uniform points in [0,1]
  Then max(y_spline) <= max(y) + 1e-14 and min(y_spline) >= min(y) - 1e-14
  And the same build with bc(1)=BC_NATURAL undershoots min(y) by more than 1e-3
      (regression guard: proves the test data actually exercises the limiter)

Scenario: Steffen limiting is third-order on smooth monotone data
  Given y = exp(x) on 33, 65, 129 and 257 uniform nodes in [0,1]
  And bc(1)=BC_STEFFEN, order(1)=3
  When evaluated at 1000 interior points
  Then the sup-norm error ratio between successive refinements is between 7.0 and 9.0
  And the same build with bc(1)=BC_NATURAL gives a ratio between 14.0 and 18.0
      (documents the deliberate loss of one order)

Scenario: quintic rejects the limiter
  Given order(1)=5 and bc(1)=BC_STEFFEN
  When construct_batch_splines_3d is called
  Then it error stops with a message naming order 3 as the only limited order

Implementation plan

Do these as two separate PRs against two separate issues; they share no code.

Issue A — compact coefficient storage (do this first; it is the one SIMPLE#423 is blocked on)

1. src/interpolate/batch_interpolate_types.f90: add to BatchSplineData3D (lines 43-58) a logical :: compact = .false. and real(dp), dimension(:,:,:,:,:,:,:), allocatable :: coeff_compact with layout (nq, 0:o1/2, 0:o2/2, 0:o3/2, n1, n2, n3) — the moment set {f, f_xx} per axis for cubic, {f, f_xx, f_xxxx} for quintic. Default .false. keeps every existing consumer byte-identical.
2. src/interpolate/batch_interpolate_3d.f90: add an optional, intent(in) :: compact argument to construct_batch_splines_3d (line 45), construct_batch_splines_3d_legacy (line 64), construct_batch_splines_3d_lines (line 207), construct_batch_splines_3d_resident (line 532), construct_batch_splines_3d_resident_device (line 543) and construct_batch_splines_3d_resident_device_impl (line 582).
3. Phase A1 (correctness): implement dense_to_compact_3d in the same module. The stored dense coefficients are Taylor coefficients about the node in the physical local coordinate — spline3d_many_point_body.inc:37-39 sets x_local = (x_norm - i)*h_step, so coeff(k,...) is f^(k)/k!. Extraction is therefore exact and trivial: f = coeff(iq,0,0,0,...), f_xx = 2*coeff(iq,2,0,0,...), f_xxyy = 4*coeff(iq,2,2,0,...), f_xxyyzz = 8*coeff(iq,2,2,2,...); for quintic add the k=4 terms with factor 4! = 24. Build dense, convert, deallocate(spl%coeff). This proves the reconstruction before touching the builders.
4. Phase A2 (the actual memory win): direct construction that never materialises the dense array. Apply the existing 1D line solve axis by axis to a working set that starts as 1 field (f) and doubles per axis: pass over x1 yields {f, f_xx}; pass over x2 applied to both yields {f, f_yy, f_xx, f_xxyy}; pass over x3 yields the 8 (cubic) / 27 (quintic) node quantities. Peak memory equals the compact array. Reuse spl_build_line_inplace (acc_spline_build_order5.f90:17) — its output work(:, line, 2) is the second Taylor coefficient, i.e. f''/2!, so the moments are already there and no new solver is needed.
5. New file src/interpolate/spline3d_compact_point_body.inc: cell reconstruction. For each axis, evaluate the two (cubic) / three (quintic) 1D basis functions at the local coordinate, then contract the 2^3 / 3^3 node quantities at the two bracketing nodes per axis — 4^3 = 64 node values for tricubic, 6^3 = 216 for triquintic. Guard with !$acc routine seq to match the existing evaluators (pattern established by 83230df, "interpolate: make batch spline evaluators acc routine seq (interpolate: make batch spline evaluators acc routine seq #283)").
6. src/interpolate/batch_interpolate_3d.f90: branch on spl%compact at the top of evaluate_batch_splines_3d_many (line 925), _many_resident (line 976), _many_der (1023), _many_der2 (1040), _der (1058), _der2 (1229), _der3 (1239). _der2_rmix (1342) and _der3 may initially error stop "not supported with compact storage" — declare that explicitly rather than silently falling back.
7. Same file: destroy_batch_splines_3d (line 887) and destroy_batch_splines_3d_device_only (line 900) must acc exit data / deallocate coeff_compact; the acc_is_present(spl%coeff) residency probes at lines 713, 892, 905 and 951 need the compact counterpart or they will silently take the host path on GPU.
8. Tests: new test/source/test_batch_interpolate_compact.f90 implementing scenarios 1-2. Register beside the existing batch tests — add_executable near test/CMakeLists.txt:187, add_test near line 276. Assert the coefficient count, so a future dense-then-convert regression is caught by CI rather than by hand.
9. Benchmark: extend bench/spline_many/ (src/bench_spline3d_many.f90) with a compact-vs-dense evaluation timing. Compact trades flops for memory traffic; on GPU it should win, on CPU it may not. Record the number in the PR — do not assume.

Issue B — Steffen-limited cubic axis

1. src/interpolate/batch_interpolate_types.f90: add integer :: bc(3) = BC_NATURAL with parameters BC_NATURAL = 0, BC_PERIODIC = 1, BC_STEFFEN = 2. Keep the existing periodic(3) logical as the authoritative wrap flag (the evaluate .inc bodies read it) and derive bc from it, or the two will disagree.

2. New file src/interpolate/spline_steffen.f90, module spline_steffen, subroutine spl_three_steffen_line(n, h, a, b, c, d) matching the signature of spl_three_reg_line (acc_spline_build_order5.f90:60): compute secant slopes s_i, the Steffen (1990) limited node slope, then the cubic Hermite coefficients converted to the same nodal Taylor layout (b, c, d) = (f', f''/2, f'''/6) that spl_three_reg_line produces — the evaluators must see an unchanged layout. Mark !$acc routine seq, no allocation, and no inner subroutines (DESIGN.md, "Trampoline Elimination").

3. Thread bc through all three builders — this is the step most likely to be done half-way:

   • construct_batch_splines_3d_legacy (64-205), dispatch sites 152/154, 172/174, 194/196;
   • construct_batch_splines_3d_lines (207-530), dispatch sites 347-497;
   • construct_batch_splines_3d_resident_device_impl (582-885) via spl_build_line_inplace (802, 827, 852), extending spl_build_line_inplace (acc_spline_build_order5.f90:17) with a bc argument.

   error stop for BC_STEFFEN with order 4 or 5 — the limiter is not defined there and a quintic "limited" spline would be a fiction. Note that the CPU default path is _legacy and the GPU default path is _resident_device_impl, so testing only one of them proves nothing about the other.

4. No evaluate-side change: the output is still nodal Taylor coefficients in the same layout.

5. Document in the module header that a BC_STEFFEN axis is C1 only and third-order, so evaluate_batch_splines_3d_der2 along that axis is piecewise-defined and discontinuous at nodes. Downstream SIMPLE#423 option (B) reads grad|B| from _der2 and must not silently use a limited axis.

6. Tests: new test/source/test_batch_spline_steffen.f90 implementing scenarios 3-5; register in test/CMakeLists.txt.

7. Mirror the same bc plumbing into batch_interpolate_1d.f90 / _2d.f90 only if a consumer needs it — otherwise leave them and say so, to keep the PR reviewable.

Both features are opt-in with defaults that preserve current behaviour, so the downstream gate (ci/downstreams, which lists SIMPLE) should stay green without a SIMPLE-side change. SIMPLE#423 depends on Issue A, not the reverse. Confirm the SIMPLE-side build integration before relying on that — it was not verified from this repository.

Acceptance criteria

• cmake --preset default && cmake --build --preset default && cd build && ctest --output-on-failure passes with zero new failures.
• test_batch_interpolate_compact passes: compact vs dense agree to 1e-12 relative for order=[3,3,3] and 1e-11 for order=[5,5,5], values and first derivatives, mixed periodic/natural.
• Compact coefficient count is exactly 8 (cubic) / 27 (quintic) doubles per node per quantity — asserted via size(spl%coeff_compact) in the test, not measured by hand.
• Peak RSS of the Phase-A2 construction path on a 120x32x120 quintic 4-quantity case is under 0.8 GB (dense path allocates 3.19 GB for coeff alone; compact is 0.40 GB).
• test_batch_spline_steffen passes: zero overshoot beyond 1e-14 on the clamped-plateau case, natural BC demonstrably failing the same bound, and a sup-norm error ratio in [7.0, 9.0] per halving on exp(x).
• test_batch_interpolate, test_batch_interpolate_many, test_batch_interpolate_rmix, test_batch_interpolate_der3, test_batch_interpolate_oracle and test_batch_eval_oracle unchanged — the default compact = .false., bc = BC_NATURAL path must be bit-identical.
• make test-nvfortran builds and the GPU-resident compact evaluate matches the CPU compact path. For Issue B, the _resident_device_impl limiter path is covered, not only _legacy.
• Benchmark number for compact-vs-dense evaluation on both CPU and GPU recorded in the PR description.

References

• M. Steffen, A simple method for monotonic interpolation in one dimension, Astronomy & Astrophysics 239, 443 (1990) — the limiter to implement.
• F. N. Fritsch and R. E. Carlson, Monotone piecewise cubic interpolation, SIAM J. Numer. Anal. 17, 238 (1980) — the monotonicity region the Steffen slope limiter lands inside.
• H. Akima, A new method of interpolation and smooth curve fitting based on local procedures, J. ACM 17, 589 (1970) — the cheaper non-monotone alternative, worth benchmarking against Steffen if the plateau bound turns out too aggressive.
• F. Lekien and J. Marsden, Tricubic interpolation in three dimensions, Int. J. Numer. Meth. Engng 63, 455 (2005) — background on the 64-coefficients-per-cell tricubic form and the flop/memory trade against a compact per-node representation. Note libneo's dense array is a nodal tensor-product Taylor expansion built from 1D line solves, not the Lekien–Marsden 64x64 linear system; the coefficient count coincides, the construction does not.
• J. Varje et al., High-performance orbit-following code ASCOT5 for Monte Carlo simulations in fusion plasmas, arXiv:1908.02482 — the code whose interp3Dcomp (compact, per-node) vs interp3Dexpl (explicit, per-cell) split is the API being requested here. Verify the exact interp3Dcomp node-quantity set against ASCOT5 source before mirroring it.
• libneo already builds ASCOT5 5.6.1 in CI (.github/workflows/main.yml:103-115, test/scripts/test_ascot5_compare.py), but that job compares the Cerfon–Freidberg analytic equilibrium solver, not splines. It is not an oracle for this work; a compact-vs-dense oracle has to be written from scratch.
• PSPLINE / EZspline is reported to use a comparable compact per-node tricubic representation for large 3D fusion field grids. Not verified here — check before citing it as precedent in the PR.

Scope and risk

Issue A is one substantial PR (new storage component, new build path, new evaluate kernel, GPU residency, tests) and should not be merged as dense-to-compact conversion only — that adds an API with none of the memory benefit. Issue B is small but touches three builders, which is where it will go wrong.

Main risk in A: the OpenACC residency probes (acc_is_present(spl%coeff) at batch_interpolate_3d.f90:713, 892, 905, 951) silently missing the new array and degrading GPU runs to host evaluation without any error.

Main risks in B: (i) implementing the limiter in only one of the three builders, so CPU and GPU disagree; (ii) a consumer picking a limited axis and then reading _der2 across it; (iii) the issue's promise of unchanged interior accuracy, which is not deliverable — a Steffen axis is C1 and third-order.

Related: itpplasma/SIMPLE#422 and #423 (the consumers) and #325 (chartmap covariant basis, merged as 5a353ba). Earlier drafts of this comment also cited libneo#204, #229 and #230; none of those numbers appear in this repository's history and their content was not verified, so they are dropped rather than described.