Use cubic interpolation for vacuum field tables - #819
Conversation
|
Codex usage limits have been reached for code reviews. Please check with the admins of this repo to increase the limits by adding credits. |
jurasic-pf
left a comment
jurasic-pf
left a comment
There was a problem hiding this comment.
Good idea! I think cubic interpolation makes a lot of sense, especially if the performance impact is barely measurable and the accuracy impact clearly is!
Given all the results you showed, I don't think there's a good reason to keep supporting linear interpolation anymore, I'm happy to just switch to the cubic default for everything. The benchmarks are convincing, and thanks for the writeup, this made the review very straightforward.
The information is tracked in the PR commits for posterity, but I don't think all of it should be merged, just the feature. I left comments in the according places.
2 participants
External-field tables use tensor cubic interpolation in R and Z. The four-node stencil shifts inward at grid edges; tables with only two or three nodes use the degree supported by their available nodes. Toroidal-plane selection and the out-of-domain error are unchanged.
The ordinary solver uses this implementation unconditionally. The linear implementation and interpolation setting have been removed, including their input/output and serialisation changes. The experiment write-up, benchmark script and results JSON are removed from the final tree. Their evidence remains in the previous PR revision.
Regression tests check tensor polynomials of degrees one, two and three, including grid edges, small tables, multiple toroidal planes and current-weighted coil groups. CTH and Solov'ev equilibrium tests check Fourier-geometry improvement under field-table refinement. The numerical method changes the free-boundary reference values, so the affected Fortran checkpoints and output fixtures are regenerated with cubic interpolation. Existing comparison tolerances are retained.
Validation of revision
e47ab40b(including current upstreamde5166b1) on Apple Silicon/macOS:Measurements retained from the earlier PR revision
At matched plasma resolution, profiles, currents and force tolerance, cubic 51x51 gives 16 to 100 times smaller physical-point B differences than linear 101x101 against a 401x401 refinement reference from identical source coils. Independent Biot-Savart/quadrature checks verify the source fields. This measures field-table error; plasma truncation and smooth-coil representation error remain separate.
Six adjacent process-isolated timing pairs per case/grid against upstream
98aaaeagive 120/120 converged solves, including asymmetric CTH and finite-pressure W7-X. The largest original-grid median slowdown is 6.1%; no solver speedup is claimed. A separate 12-configuration, three-profile QUASR matrix has the same five of 36 converged problems with both methods at the supplied budgets. An additional W7-X check at measured 4.24% total beta converges in all four method/grid combinations; cubic 61 has about 59 times smaller physical-point B difference than linear 61 against the original 121-table cubic reference.These are historical measurements from the earlier revision. The final cubic weights on tables with at least four nodes retain its closed-form expressions.
Reference regeneration and exact Fortran source patch
The independent reference generator is educational_VMEC at ce663b2 with the patch below. It also includes the asymmetric debug-output slices from educational_VMEC#30. It changes
becoilto tensor Lagrange interpolation and restores the VMEC 8.52lforbalcalculation for the existing Solov'ev force-balance test. This is a deliberately modified reference solver, not unmodified VMEC 8.52.To validate the restoration, a build with the original bilinear
becoilreproduces the existinglforbal=truegolden's R/Z/lambda/B Fourier coefficients to a worst normalised difference of 4.8e-11. The cubic build then supplies the replacement golden.Build with CMake/Ninja, gfortran, NetCDF C/Fortran and FFTW. The local compiler flags were
-ffree-line-length-none -fallow-argument-mismatchplus the NetCDF-Fortran include directory. Runxvmec input.cth_like_free_bdyandxvmec input.solovev_free_bdyusing the inputs fromsrc/vmecpp/cpp/vmecpp_large_cpp_tests/test_dataand the unchanged mgrid files fromsrc/vmecpp/cpp/vmecpp/test_data. Those inputs contain the checkpoint dump flags. Replace matching JSON checkpoints when the maximum normalised numerical change exceeds 1e-12; retain unchanged checkpoints. The final multigrid filenames follow the new Fortran iteration numbers: CTH 496, Solov'ev 417 and 808.The large-test CTH wout and the core Solov'ev wout are Fortran outputs. For
wout_solovev_free_bdy_T.nc, use the coreinput.solovev_free_bdywithLFORBAL = T. The core CTH wout is a VMEC++ serialisation fixture, as documented in #781; regenerate it with the final VMEC++ implementation, one thread and the unchanged CTH input/mgrid, retaining the full VMEC++ output schema.The three asymmetric
vac1n_analyt/vac1n_fourireferences introduced by #738 are regenerated from the corecth_like_free_bdy_asym.json, converted withvmecpp._util.vmecpp_json_to_indata, and the unchangedmgrid_cth_like_asym.nc. Setiter2_to_dump = 53, 54and enabledump_vac1n_analyt,dump_vac1n_fourpanddump_vac1n_fouri. The geometricfourpfixture is unchanged. The bilinear reference build reproduces the original asymmetric bvec values within 3e-14; the rebuilt cubic reference passes both new asymmetric C++ targets at their existing tolerances.Run each reference case in a fresh output directory: educational_VMEC stops when a checkpoint JSON already exists, even though its process exit code can be zero. Verify
EXECUTION TERMINATED NORMALLYin the run log. Use a clean build directory when changing the reference interpolation method.The following patch applies to the reference-generator commit linked above: