Utilities - BUG - CubicFit end derivative loses an offset-sized rounding error on edge-packed grids

Open
  • Is something wrong?
#559 0 comments 0 reactions 0 assignees View on GitHub

Maintainers usually reply within 3 days

Nobody has claimed this yet.

Assessment

Difficulty
4/5
Estimated time
3-5 days
Newbie friendliness
49/100
Issue type
Bug
Clarity
Mostly clear
Activity status
Active
Tech stack
julia
Domain
performance

Research direction

Start with FastInterpolations.jl's src/core/polyfit_kernels.jl:55 and :612-622, then inspect GPEC.jl's edge-packed grid construction in DirectEquilibrium.jl:478-500. Reproduce the end-derivative error using an offset profile on an edge-packed grid and compare against the analytic derivative. If confirmed, the issue proposes fixing it upstream in FastInterpolations.jl and updating this repository's compatibility bound.

Written by the indexing model from the issue text.

Description

GPEC.jl builds its equilibrium splines with cubic_interp, whose default boundary condition CubicFit() sets the end derivative from a 4-point polynomial fit. FastInterpolations.jl 0.4.20 evaluates this derivative as muladd(c1, f1, … c4*f4) (src/core/polyfit_kernels.jl:55, called from :612-622). The weights sum to zero, so a constant offset in f cancels only in exact arithmetic. In floating point the error is about eps*|f|/h and is not reduced, so it grows on grids whose end spacing h is small compared with the end value.

The Fortran code had exactly this defect in its extrap spline end derivative. It made DCON/GPEC results wrong for mpsi >= 320 with grid_type='ldp' and at large aspect ratio (GPEC#306, fixed in GPEC#308). The fix there forms the derivative from differences to the end value, sum c_i (f_i - f_1), which is algebraically identical and removes the offset before rounding. The edge-packed Julia grids (log_asymptotic, core_edge_packed, DirectEquilibrium.jl:478-500) create the small h where this matters.

Status: verified in the source; the effect on a GPEC.jl result is not yet measured. Proposed route: (1) reproduce with a profile f = C + g(psi) that has a large offset C, on an edge-packed grid, and compare the end derivative with g' analytically; (2) if confirmed, repair it upstream in FastInterpolations.jl and bump the compat bound here, as for the PeriodicBC fix (#531).

Chris&AI

Dominant language
Julia
Stars
11
Forks
7
Avg merge
3d 22h
Merged PRs (30d)
70

Getting set up

  • No Dockerfile or Docker Compose file
  • Has a pull request template
  • No contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from OpenFUSIONToolkit/GPEC

All issues in OpenFUSIONToolkit/GPEC

Similar issues

More Julia issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.