celeritas-project
Implement spline interpolation
Spline interpolation (or some other nonlinear interpolation) is used in Geant4 by :
- Energy loss rate calculation (see
calc_energy_loss) - LivermorePE cross sections above K-shell energy (see LivermoreElement)
- Cross section (mean free path) for multiple scattering models (Urban and WenzelVI)
Research why it's necessary (the energy loss one provides an explanation of the consequences of linear interpolation, although it seems to me an alternate interpolation scheme such as log-log might have the same effect with less overheard) and compare spline interpolation with alternatives.
@amandalund and @stognini : you ran some instances of TestEm3 with and without spline interpolation, and with refined energy grids. Can you post/reference any of the results/analysis when you get a chance?
Here is the test I currently have. Celeritas is showing slightly larger differences. Plot details:
- G4 v10.7.2; binning was not changed
- 100k primaries, 1mm threshold cuts for all samples
- Relative differences are
(spline_plot – X) / spline_plot
I had run TestEm3 with Celeritas with refined energy grids (1e7 primaries, 1mm production cuts):
Default is 7 bins/decade, and the error is relative to the 112 bins/decade result.
Just ran a short test on Summit. The input, cpu run script, and results are in benchmarks/testem3-demo-loop/g4app-summit/, but here is the gist of it:
| Spline | Avg. wall time [s] |
|---|---|
| ON | 176.2 |
| OFF | 170.5 |
Time penalty is roughly 3%.
Nice! I'll do a run with increased data points for comparison to see what the penalty is there.
I've added an interface to our input to allow adjustment of the number of data points: see https://github.com/celeritas-project/celeritas/pull/414 .
Let's wait until refactoring the physics to use standard 'builders' before doing this.
@amandalund @lebuller Moving https://github.com/celeritas-project/celeritas/pull/1562#discussion_r1927796966 here. I thought spline interpolation is basically a local polynomial interpolation, which is what #1444 implements? So perhaps y'all could make sure we're on the same page and able to reproduce Geant4's interpolation method.
I think they are very similar, but cubic splines have the constraints that the first and second derivatives be continuous everywhere? A downside is this would require storing the second derivatives as well, and what we have might be close enough: it looks like the errors due to the differences in the interpolation method In #1544 at least are pretty small. Hacking together something equivalent to what Geant4 does (which looks like cubic spline interpolation with not-a-knot boundary conditions) and using it to build the range table, I'm able to exactly reproduce the range table values from Geant4.
For context, some nice examples/theory at https://blog.timodenk.com/cubic-spline-interpolation/index.html . I guess it would behoove us (@lebuller ) to do a little convergence study of the range (which is the most important integral quantity) with the different methods: increasing geant4 xs points, then we compare the imported range (i.e. using their not-a-knot interpolation) to the calculated range (using your calculator here with the cubic spline eloss?)