Skip to content

Preserve constant offsets in spline and cspline endpoint derivatives - #308

Merged
matt-pharr merged 4 commits into
PrincetonUniversity:developfrom
krystophny:fix/spline-extrap-endpoint-cancellation
Oct 16, 2026
Merged

matt-pharr merged 4 commits into
PrincetonUniversity:developfrom
krystophny:fix/spline-extrap-endpoint-cancellation

Conversation

@krystophny

Copy link
Copy Markdown
Collaborator

Fixes #306.

spline_fit_ahg and cspline_fit_ahg (endmode "extrap") computed the end derivative as SUM(c*f) with
c(1) = -sum(c(2:4)). The terms are O(f/h) and cancel, so a constant or large-offset profile on a grid clustered at
the right end acquires an error of order eps*|f|/h (4.8e-2 for h ~ 1e-5, |f| ~ 1). The change evaluates the same
four-point derivative from differences to the end value, in equil/spline.f and equil/cspline.f. No other
behaviour changes.

Test: regression/test/test_spline_endpoint_polynomials.py fits constant and cubic profiles (and a complex constant
profile with offset 300) on the 'ldp' sine-squared grid with 512 and 1024 intervals and checks the end derivatives,
the end values and the integral through spline_int. The test compiles only io, local, utils, spline and
cspline from equil/. Before the change it fails at 512 and 1024 intervals for every profile (largest derivative
errors: 0.11 and 12.4 for the cubic, 0.89 for the constant, 586 for the complex constant); after the change it passes.

End to end (lar deck of the issue, A = 40, 3/1 drive): DCON plasma energy and GPEC response are the same for
mpsi = 256, 320, 384, 512 (energy -6.434e-2, difference below 1e-3 relative) and for ma = 512, 1024; for
A = 70..1000 the total free-boundary energy is positive and scales as 1/A.

Chris&AI

The 'extrap' end derivatives were evaluated as SUM(c*f) with c(1) = -sum(c(2:4)).
For a constant profile, or any profile with a large offset, on a grid clustered
at an endpoint the four terms cancel catastrophically (error ~ eps*|f|/h^3).
Evaluate the same four-point derivative from differences to the endpoint value.

@logan-nc logan-nc left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds good to me and I know @matt-pharr has issues with the ldp grid packing in the past. I haven't independently checked this though.

@logan-nc

Copy link
Copy Markdown
Contributor

@matt-pharr I wanted to ping you on this since I know you've had ldp spline packing issues in the past. You are assignee just to make sure you see it. You can hit merge or do a more thorough review - challenge by choice. I was convinced by the included test but didn't independently check.

@matt-pharr

Copy link
Copy Markdown
Collaborator

@logan-nc it makes me very uncomfortable that we are changing the splines... I will need to look into this very carefully before merging

@matt-pharr

Copy link
Copy Markdown
Collaborator

@krystophny FYI a bunch of your recent PRs have been into master, they all should be into develop, and you should be using develop in general as that is the most up to date version. Master is for releases.

@matt-pharr
matt-pharr changed the base branch from master to develop October 15, 2026 10:02
@krystophny

Copy link
Copy Markdown
Collaborator Author

@krystophny FYI a bunch of your recent PRs have been into master, they all should be into develop, and you should be using develop in general as that is the most up to date version. Master is for releases.

Understood. I will make Claude understand that too.

@matt-pharr

matt-pharr commented Oct 16, 2026 •

Copy link
Copy Markdown
Collaborator

@logan-nc I had claude run some tests on the errors with the spline packing, and it seems like this hearkens way back to the classic spline investigations by @zrwang2015... #68 #55. Here are some test representations of different functions with right-hand-side packing:

case mpsi original AHG @krystophny patched AHG @zrwang2015 classic
exp(2x)cos(5x) 256 7.0e-4 1.0e-7 1.0e-7
exp(2x)cos(5x) 512 4.4e-2 6.4e-9 6.4e-9
exp(2x)cos(5x) 1024 3.6 4.0e-10 4.0e-10
f = 40 + x²/2 (F-like offset) 512 0.86 6.8e-10 1.3e-9
f = 1 1024 1.7 0 0
complex 300−2i 1024 521 0 —
exact cubic 1024 24 2.2e-10 5.0e-10

It also found the reason that zhirui's splines were producing different/worse results in some cases, that being that it uses real types where it should use complex types in cspline_fit_classic and therefore throws away complex derivatives. So I am going to work that fix into this PR. I also checked out chris' changes and they do seem correct. I suggest merging this after I clean up the extraneous python tests.

Funny enough this change would have saved me from the year of despair at the RDCON 'infinite surface shielding' issue :')

@matt-pharr matt-pharr left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed the tests and added the type fix for the classic splines. This is merge-able now.

@matt-pharr
matt-pharr merged commit 030a102 into PrincetonUniversity:develop Oct 16, 2026
6 checks passed
@logan-nc

Copy link
Copy Markdown
Contributor

Thanks @krystophny for this deep fix and @matt-pharr for the detailed checking / tracking. This is definitely great to have, and hopefully no more researchers lose a year to faulty spline!

The classic spline fix is great to have too. But those ate definitely slower, so this one is a value add beyond that 😀

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

3 participants