Repository navigation
EQUIL, RDCON, PENTRC - FEATURE - Monotone (pchip) fits for the tabulated Zeff and kinetic profiles - #300
Conversation
…bulated data without derivatives Fritsch-Carlson interior slopes with scipy PchipInterpolator's end slopes, so Python front ends (bouquet, TPS) can reproduce GPEC's interpolant exactly. Unit test checks values and derivatives against scipy. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…profile table The tabulated Zeff (RDCON MRE terms) and the PENTRC kinetic input table carry pedestal-scale gradients; a cubic spline rings there, pchip does not. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
pchip fits any table of 2 or more points, so the 4-point minimum of the "extrap" spline goes. A scalar Zeff now becomes a flat profile on 2 nodes, the same table that tearing-physics-suite writes for a scalar Zeff (psi_N_Zeff = 0, 1). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
@StuartBenjamin this touches the splines, so I am going to need a bit to review this because that is the well from which we all drink :) |
|
This actually mostly just creates a new class of spline so it might be easy, and its main implications are for PENTRC which @logan-nc might have some thoughts on. |
|
Skimming - I like the concept of this type of spline for the input profiles, and I like the minimal highly-independent implementation. @StuartBenjamin, can you post some plots showing these ringings you describe in the top description and the fix? These would be very helpful now and for future reference. |
Conflicts: equil/spline.f (develop's spline_fit_hermite stays subprogram 24; spline_fit_pchip becomes 25) and rdcon/dcon.f (keep >= 2 Zeff knots with the 2-point default table, plus develop's strictly-increasing check). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
In short the pchip doesn't make a huge difference unless you've got a sparse or bumpy input 
The PENTRC thing that the robot identified in the PR report is a separate issue (see #321 fresh and ready to be investigated): I wrote this PR because I wanted a bulltproof spline that I didn't have to worry about, it's marginal gain if you have well behaved profiles already |
…rozen for v1.0) The test passed locally against scipy's PchipInterpolator before removal; its source remains in the PR #300 description and scripts zip. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
logan-nc
left a comment
logan-nc
left a comment
There was a problem hiding this comment.
Approving. I went through spline_fit_pchip line by line against scipy's _find_derivatives and _edge_case. The interior weights, the zero slope at extrema, both end-slope branches (the d0 == 0 check reproduces numpy's sign(0) = 0) and the 2-point secant all match. Before removing the unit test, I built it locally and got PASS, and I recomputed its reference values with scipy.
I pushed f16a602, which removes regression/spline_tests/. Standalone tests and new test plumbing are frozen until the Fortran v1.0 release (see d4160a9 and 89d53a7). The test source is still in the PR description and the scripts zip.
A note for the record: in read_kin, pchip fits only the intermediate table spline. kin is still resampled out to ψ_N = 1, which is past the table end, and refit with a cubic. So the "no overshoot" guarantee does not reach what PENTRC uses at the edge, and the ±21–25 % δW_k shift in the CI Solovev kinetic case is the edge artifact in #321, not pchip itself. #321 should be fixed separately.
🤖 Generated with Claude Code
Builds on #296. Fits the tabulated Zeff profile (RDCON MRE terms) and the PENTRC kinetic profile table (also read by DCON with
kin_flag) with a monotone cubic (pchip) instead of a cubic spline. These tables have pedestal-scale gradients and no derivatives. A cubic spline overshoots and rings there; pchip cannot overshoot between points.Changes
spline_fit_pchip(equil/spline.f): monotone cubic Hermite fit with the same slopes as scipy'sPchipInterpolator, for 2 or more points.rdcon/mercier.f) and inread_kin(pentrc/inputs.f90).regression/spline_tests/test_spline_pchip.f90(manual, not in CI): matches scipy to 1e-7 (values) and 1e-6 (slopes); 2-point tables are exact.Effect on results
kin_flag)mufrac1.9 %, Δ′ identicalThe PENTRC changes are an edge artifact of the example
kinetic.txt, which ends at ψ_N = 0.998 with T = 0.59 eV. PENTRC samples the table fit on a regular grid up to ψ_N = 1. There T is 1.4e-6 eV with the cubic spline and 1e-4 eV with pchip, so logΛ < 0 and ν is huge and negative at that node in both cases. PENTRC's regular-grid spline of ν (still a cubic spline) rings from that node far inward. pchip makes the edge value ~1000× smaller, so the ringing shrinks (ν_i at ψ_N = 0.87: 1.4e9 before, 3.9e6 now, against 1.5e6 in the core), but ν stays unphysical near the edge either way. n, T and ω_E change by ~1e-6.Scripts: pchip_profiles_scripts.zip
run_split_tests.sh OUTDIRruns the cases above (edit the build paths at the top), andcmpnc.pycompares the netCDF files.Dependencies
Contains #296. Conflicts with #297 only where both append code (
equil/spline.f: keep both,spline_fit_pchipbecomes subprogram 25; the spline-test makefile and .gitignore: build both tests).🤖 Generated with Claude Code