Skip to content

Name first generated kernel: splined field evaluation (#515) #516

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
krystophny wants to merge 1 commit into
base: main
Choose a base branch
from
+339 −0
Open

Name first generated kernel: splined field evaluation (#515) #516

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
133 changes: 133 additions & 0 deletions DOC/first-generated-kernel.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,133 @@
# First Generated Kernel: Splined Field Evaluation

This document records the decision that closes the gap between the
lazy-fortran generation programme and a production SIMPLE physics path
(GitHub issue #515). The lazy-fortran stack generates numerical kernels from
symbolic definitions and emits them to CPU Fortran, OpenMP offload, OpenACC and
native CUDA (tracker: lazy-fortran/fortgen#1). Every issue in that stack ends
inside the toolchain itself; none of them touches code that publishes a physics
number. This page names the first SIMPLE kernel that a generated implementation
must replace, explains why, and fixes the measurement instrument that decides
whether the generated kernel is worth shipping.

The decisive sentence this document sets out to make:

> For this kernel, the arithmetic the mathematics requires is `N_sym`, the
> generated code emits `N_emit`, the compiler produces `N_machine`, the
> hardware delivers `T`, and the accuracy is `X` ulp against a high-precision
> reference.

---

## 1. The choice: splined field evaluation

**The first generated consumer is `splined_evaluate` / `splined_evaluate_with_der`
in `src/field/field_splined.f90`.**

These two subroutines are the leaf of the per-particle orbit-tracing path. Every
guiding-centre step that uses a splined field calls them, they are straight-line
and arithmetic-heavy (no control flow beyond the batch-spline basis sum), and the
derivative variant `splined_evaluate_with_der` has exactly the shape symbolic
automatic differentiation produces well. This makes it the natural first target
for the generation programme.

### Why not the alternatives

1. **Canonical-coordinate transformation kernels** (`field_can_*` family):
larger, with substantially more control flow and coordinate-system dispatch.
A poorer first target: the symbolic source is bigger, the emitted code is
harder to verify, and the win per unit of toolchain effort is smaller.
2. **RK stage evaluation**: already generated inside the toolchain (fortnum),
but SIMPLE does not consume the generated version, and the RK stages call the
field anyway, so the field kernel is upstream of them.

### Evidence that this is where the time goes

- The spline path is the hottest evaluation point in the code: it replaces the
Biot-Savart / VMEC field evaluation on every orbit step.
- libneo#408's 2.16x spline evaluation improvement is direct evidence that this
path dominates wall time; halving its data (fortnum#84, mixed precision) is
flagged as the largest available win, consistent with a memory-bound
straight-line kernel.
- The batch-spline optimization already merged into SIMPLE (v1.4.0, 2.86x
speedup) measured the same routine.

The choice is a decision recorded here, not an inference from file size: the
baseline instrument below measures `splined_evaluate` on a real NCSX workload,
so the claim "this is the hot path" is testable in-repo.

---

## 2. What "first generated consumer" means

A generated implementation of `splined_evaluate_with_der` must be drop-in
compatible: same argument list, same `BatchSplineData3D` coefficients, same
outputs, evaluated at the same points, and interchangeable with the hand-written
subroutine behind the `splined_field_t` type. The generated kernel competes with
the hand-written one on the four emitted targets (CPU Fortran, OpenMP offload,
OpenACC, native CUDA), and the acceptance gate is:

- **Golden records unchanged.** The numerical gate is the existing one; a
generated kernel that changes any golden-record output is rejected regardless
of speed.
- The kernel is verified by readback against its symbolic source
(lazy-fortran/fortnum#70), checked against the hand-written implementation
numerically, and measured with the instrument in this document.

---

## 3. The measurement instrument

`test/tests/test_splined_kernel_bench.f90` measures the hand-written kernel on a
realistic NCSX workload (coils field + VMEC reference coordinates). It reports
the two numbers SIMPLE can measure itself:

- **`T`**: wall-clock time per evaluation of `splined_evaluate` (value only) and
`splined_evaluate_with_der` (value + derivative), best-of-N over repeated
loops of `n_timing` evaluations.
- **`X` ulp**: maximum ulp deviation of the splined-field `Bmod` from the direct
Biot-Savart coils field at the same physical point, across a sweep of points.
The direct coils field is the high-precision reference for this kernel. The
reported figure is dominated by spline interpolation error (the spline is an
approximation of the reference field), which is exactly the accuracy the
kernel delivers as a physics producer.

The other two counts come from the toolchain, not from SIMPLE:

- **`N_sym`**: the arithmetic the mathematics requires, counted in the symbolic
definition in fortsym.
- **`N_emit`**: the arithmetic the generated code emits, counted by fortsym.
- **`N_machine`**: the machine instructions the compiler produces, counted from
the generated code (e.g. assembly instruction count).

When a generated kernel lands, the instrument is rerun unchanged and the
generated numbers are inserted into the statement. The statement is the
deliverable: it is the only way to tell whether the generated kernel is faster,
slower, or as accurate as the hand-written one, and that finding belongs in
SIMPLE where the consequence lands.

### Baseline (this repository, gfortran 14, one core)

Measured by `test_splined_kernel_bench.x` on the NCSX coils workload:

```
T (value only) = ~1.1e-6 s/eval
T (value + derivative) = ~3.2e-6 s/eval
X ulp (max, Bmod) = ~2.7e9 ulp vs direct field (~6e-7 relative)
```

The ulp figure is spline-interpolation accuracy at the default 62x63x64 grid,
not a rounding-error bound. It is the accuracy the kernel delivers today; a
generated kernel must match it (or the golden-record gate rejects it).

---

## 4. Related

- lazy-fortran/fortgen#1 -- tracker for the generation and measurement programme
- lazy-fortran/fortnum#78 -- end-to-end objective `benchmark-orbit-proxima`
- lazy-fortran/fortnum#70 -- readback verification against the symbolic source
- lazy-fortran/fortnum#84 -- mixed precision; the spline path is memory-bound
- lazy-fortran/fortnum#80 -- attribution of the 1.70x CUDA-over-OpenACC gap
- `src/field/field_splined.f90` -- the kernel this document names
- `test/tests/test_splined_kernel_bench.f90` -- the measurement instrument
10 changes: 10 additions & 0 deletions test/tests/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -745,6 +745,16 @@ add_test(NAME test_bder_from_splines
set_tests_properties(test_bder_from_splines PROPERTIES
LABELS "integration")

add_executable(test_splined_kernel_bench.x test_splined_kernel_bench.f90)
target_link_libraries(test_splined_kernel_bench.x simple)
add_dependencies(test_splined_kernel_bench.x ncsx_test_data)
add_test(NAME test_splined_kernel_bench
COMMAND test_splined_kernel_bench.x
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
set_tests_properties(test_splined_kernel_bench PROPERTIES
LABELS "unit;performance"
TIMEOUT 120)

add_executable(test_refcoords_file_detection.x test_refcoords_file_detection.f90)
target_link_libraries(test_refcoords_file_detection.x simple)
add_dependencies(test_refcoords_file_detection.x chartmap_test_data)
Expand Down
196 changes: 196 additions & 0 deletions test/tests/test_splined_kernel_bench.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,196 @@
program test_splined_kernel_bench
!> Baseline measurement instrument for the first generated consumer
!> (splined field evaluation, see DOC/first-generated-kernel.md).
!>
!> This program measures the hand-written splined field kernel on a
!> realistic NCSX workload and reports the two numbers that SIMPLE itself
!> can measure for the fortsym statement:
!>
!> - T: wall-clock time per evaluation of splined_evaluate and
!> splined_evaluate_with_der (the hottest path in the code), and
!> - X ulp: accuracy of the splined field path against the direct
!> Biot-Savart coils field, which is the high-precision reference
!> for this kernel.
!>
!> The arithmetic-count numbers N_sym (mathematics) and N_emit (generated
!> code) come from the lazy-fortran/fortsym generation toolchain, and
!> N_machine from the compiler; they are filled in when a generated kernel
!> is substituted for the hand-written one. This program is the
!> production-side gate the generated kernel must beat.
!>
!> Golden records are not touched: the measurement is informational and
!> only sanity-checks that the spline path produces finite values within
!> the spline-interpolation accuracy of the reference field.

use, intrinsic :: iso_fortran_env, only: dp => real64
use simple, only: init_vmec
use util, only: twopi
use new_vmec_stuff_mod, only: nper
use field_coils, only: coils_field_t, create_coils_field
use field_splined, only: splined_field_t, create_splined_field
use reference_coordinates, only: init_reference_coordinates, ref_coords

implicit none

type(coils_field_t) :: raw_coils
type(splined_field_t) :: splined

real(dp) :: dummy
real(dp) :: phi_period

integer, parameter :: n_timing = 20000 ! evaluations for the timing loop
integer, parameter :: n_rep = 5 ! repeats, best (min) time kept
integer, parameter :: n_acc = 400 ! points for the accuracy sweep

real(dp) :: t_eval, t_eval_der
real(dp) :: max_ulp

call init_vmec('wout_ncsx.nc', 5, 5, 5, dummy)
call init_reference_coordinates('wout_ncsx.nc')

call create_coils_field('coils.simple', raw_coils)
call create_splined_field(raw_coils, ref_coords, splined)

phi_period = twopi/real(nper, dp)

! --- T: end-to-end timing of the hand-written splined field kernel ---
call time_kernel(splined, phi_period, n_timing, n_rep, t_eval, t_eval_der)

! --- X ulp: accuracy of the spline path vs the direct coils field ---
call measure_accuracy(raw_coils, splined, phi_period, n_acc, max_ulp)

print *, ''
print *, '=================================================================='
print *, 'Splined field kernel baseline instrument (DOC/first-generated-kernel.md)'
print *, '------------------------------------------------------------------'
print *, 'kernel: splined_evaluate / splined_evaluate_with_der'
print *, 'reference: direct Biot-Savart coils field (wout_ncsx.nc)'
print *, '------------------------------------------------------------------'
print '(a,es12.4,a)', ' T (value only) = ', t_eval, ' s/eval'
print '(a,es12.4,a)', ' T (value + derivative) = ', t_eval_der, ' s/eval'
print '(a,es12.4,a)', ' X ulp (max, Bmod) = ', max_ulp, ' ulp vs direct field'
print *, '------------------------------------------------------------------'
print *, 'statement: for this kernel the arithmetic the mathematics'
print *, 'requires is N_sym (fortsym), the generated code emits N_emit'
print *, '(fortsym), the compiler produces N_machine, the hardware'
print *, 'delivers T above, and the accuracy is X ulp above against a'
print *, 'high-precision reference.'
print *, '=================================================================='

! Sanity gate: finite, positive timing and finite spline accuracy. The
! spline is an approximation of the reference field, so the ulp figure is
! the spline-interpolation accuracy, not a rounding-error bound; we only
! require it to be finite and reasonable (< 1e-3 relative).
if (.not. (t_eval > 0.0_dp .and. t_eval_der > 0.0_dp)) then
print *, 'FAILED: non-positive kernel timing'
error stop 1
end if
if (.not. (max_ulp >= 0.0_dp .and. max_ulp < 1.0e13_dp)) then
print *, 'FAILED: spline accuracy out of expected range'
error stop 1
end if

print *, 'PASSED: splined kernel baseline instrument'
end program test_splined_kernel_bench


subroutine time_kernel(f, phi_period, n_timing, n_rep, t_eval, t_eval_der)
!> Best-of-n wall-clock time per evaluation of splined_evaluate and
!> splined_evaluate_with_der. Outputs are accumulated into a checksum so
!> the compiler cannot elide the kernel calls.
use, intrinsic :: iso_fortran_env, only: dp => real64
use field_splined, only: splined_field_t

type(splined_field_t), intent(in) :: f
real(dp), intent(in) :: phi_period
integer, intent(in) :: n_timing, n_rep
real(dp), intent(out) :: t_eval, t_eval_der

real(dp) :: x(3), Acov(3), hcov(3), Bmod
real(dp) :: dAcov(3, 3), dhcov(3, 3), dBmod(3)
real(dp) :: checksum, dt
integer(selected_int_kind(18)) :: t0, t1, cr
integer :: i, rep

t_eval = huge(1.0_dp)
t_eval_der = huge(1.0_dp)

! Warm up
x = [0.5_dp, 1.0_dp, 1.0_dp]
call f%evaluate(x, Acov, hcov, Bmod)
call f%evaluate_with_der(x, Acov, hcov, Bmod, dAcov, dhcov, dBmod)

do rep = 1, n_rep
checksum = 0.0_dp
call system_clock(count=t0)
do i = 1, n_timing
x(1) = 0.02_dp + 0.96_dp*real(mod(i*7, 1000), dp)/1000.0_dp
x(2) = phi_period*real(mod(i*13, 1000), dp)/1000.0_dp
x(3) = phi_period*real(mod(i*29, 1000), dp)/1000.0_dp
call f%evaluate(x, Acov, hcov, Bmod)
checksum = checksum + Bmod
end do
call system_clock(count=t1, count_rate=cr)
dt = real(t1 - t0, dp)/real(max(cr,1), dp)
t_eval = min(t_eval, dt/real(n_timing, dp))

checksum = 0.0_dp
call system_clock(count=t0)
do i = 1, n_timing
x(1) = 0.02_dp + 0.96_dp*real(mod(i*7, 1000), dp)/1000.0_dp
x(2) = phi_period*real(mod(i*13, 1000), dp)/1000.0_dp
x(3) = phi_period*real(mod(i*29, 1000), dp)/1000.0_dp
call f%evaluate_with_der(x, Acov, hcov, Bmod, dAcov, dhcov, dBmod)
checksum = checksum + Bmod
end do
call system_clock(count=t1, count_rate=cr)
dt = real(t1 - t0, dp)/real(max(cr,1), dp)
t_eval_der = min(t_eval_der, dt/real(n_timing, dp))
end do

! Keep the checksum observable to the optimizer.
if (checksum == huge(1.0_dp)) print *, checksum

end subroutine time_kernel


subroutine measure_accuracy(raw, splined, phi_period, n_acc, max_ulp)
!> Max ulp deviation of splined Bmod from the direct coils Bmod at the
!> same physical point. Splined coordinates are (r=sqrt(s), theta, zeta);
!> the physical point is obtained by s = r^2 and ref_coords%evaluate_cart.
use, intrinsic :: iso_fortran_env, only: dp => real64
use field_coils, only: coils_field_t
use field_splined, only: splined_field_t
use reference_coordinates, only: ref_coords
use util, only: twopi

type(coils_field_t), intent(in) :: raw
type(splined_field_t), intent(in) :: splined
real(dp), intent(in) :: phi_period
integer, intent(in) :: n_acc
real(dp), intent(out) :: max_ulp

real(dp) :: x(3), u(3), xcart(3)
real(dp) :: Acov(3), hcov(3), Bmod
real(dp) :: Acov_ref(3), hcov_ref(3), Bmod_ref
real(dp) :: ulp
integer :: i

max_ulp = 0.0_dp
do i = 1, n_acc
x(1) = 0.05_dp + 0.9_dp*real(mod(13*i, 997), dp)/997.0_dp
x(2) = 0.1_dp + (twopi - 0.2_dp)*real(mod(29*i, 991), dp)/991.0_dp
x(3) = 0.1_dp + (phi_period - 0.2_dp)*real(mod(31*i, 983), dp)/983.0_dp

call splined%evaluate(x, Acov, hcov, Bmod)

u(1) = x(1)**2
u(2) = x(2)
u(3) = x(3)
call ref_coords%evaluate_cart(u, xcart)
call raw%evaluate(xcart, Acov_ref, hcov_ref, Bmod_ref)

ulp = abs(Bmod - Bmod_ref)/max(spacing(Bmod_ref), tiny(1.0_dp))
max_ulp = max(max_ulp, ulp)
end do
end subroutine measure_accuracy
Loading