Numeric grids¶

Structured data is used in cross section evaluation and sampling routines.

Grid access¶

Grids classes act like vectors, allowing direct access to node-centered data.

class UniformGrid¶

Interact with a uniform grid of increasing values.

This simple class is used by physics vectors and classes that need to do lookups on a uniform grid.

Searching for a grid index with find requires the input value to be in range, because out-of-bounds values often need special treatment (e.g., clipping to the boundary values rather than interpolating). It’s easier to check for those exceptional cases outside of the grid view.

Public Functions

inline explicit UniformGrid(UniformGridData const &data)¶: Construct with data.

inline size_type size() const¶: Number of grid points.

inline value_type front() const¶: Minimum/first value.

inline value_type back() const¶: Maximum/last value.

inline value_type operator[](size_type i) const¶: Get the value at the given grid point.

inline size_type find(value_type value) const¶: Find the value bin such that data[result] <= value < data[result + 1].

template<class T> class NonuniformGrid¶

Interact with a nonuniform grid of increasing values.

This should have the same interface (aside from constructor) as UniformGrid. The grid is closed on the lower bound and open on the upper, and the same treatment is used for coincident grid points: the higher of the points will be chosen.

These grids are often used with a helper function for safely finding points that can be interpolated inside the grid’s boundaries.

template<class Grid> inline auto celeritas::find_interp(Grid const &grid, typename Grid::value_type value)¶

Find the index of the value and its fraction between neighboring points.

The grid class should have a floating point value and must have methods find, front, back, and operator[] .

The value must be bounded by the grid and less than the final value. The result will always have an index such that its neighbor to the right is a valid point on the grid, and the fraction between neghbors may be zero (in the case where the value is exactly on a grid point) but is always less than one. If the requested point is exactly on a coincident grid point, the lower point and a fraction of zero will result.

Local interpolation¶

At the core of between-point evaluation is interpolation, which can be done as linear, logarithmic, semi-log, or spline.

template<Interp XI = Interp::linear, Interp YI = Interp::linear, typename T = ::celeritas::real_type> class Interpolator¶

Interpolate, with either linear or log in x and y.

The inputs are given as two (x,y) pairs. The same two pairs can be used to interpolate successively with this functor.

Interpolation on the transformed coordinates is the solution for \( y \) in

\[ \frac{f_y(y) - f_y(y_l)}{f_y(y_r) - f_y(y_l)} = \frac{f_x(x) - f_x(x_l)}{f_x(x_r) - f_x(x_l)} \]

where \( f_d(v) = v \) for linear interpolation on the \( d \) axis and \( f_d(v) = \log v \) for log interpolation.

Solving for \(y\), the interpolation can be rewritten to minimize transcendatal operations and efficiently perform multiple interpolations over the same point range:

\[ y = f^{-1}_y \left( f_y(y_l) + \frac{ p_y(n_y(y_l), y_r) }{ p_x(n_x(x_l), x_r)} \times p_x(n_x(x_l), x) \right) \]

where transformed addition is \(p_y(y_1, y_2) \equiv f_y(y_1) + f_y(y_2)\) , transformed negation is \(n_y(y_1) \equiv f^{-1}_y( - f_y(y_1) )\) and \( f_y(y) = y \) for linear interpolation in y and \( f_y(y) = \log y \) for log interpolation in y.

Instantiating the interpolator precalculates the transformed intercept and slope terms, as well as the negated x-left term. At each evaluation of the instantiated Interpolator, only the inverse-transform and add-transformed operation need be applied.

Template Parameters:

XI – Transform to apply to X coordinate
YI – Transform to apply to Y coordinate
T – Floating point type

template<typename T = ::celeritas::real_type> class SplineInterpolator¶

Interpolate using a cubic spline.

Given a set of \( n \) data points \( (x_i, y_i) \) such that \( x_0 < x_1 < \dots < x_{n - 1} \), a cubic spline \( S(x) \) interpolating on the points is a piecewise polynomial function consisting of \( n - 1 \) cubic polynomials \( S_i \) defined on \( [x_i, x_{i + 1}] \). The \( S_i \) are joined at \( x_i \) such that both the first and second derivatives, \( S'_i \) and \( S''_i \), are continuous.

The \( i^{\text{th}} \) piecewise polynomial \( S_i \) is given by:

\[ S_i(x) = a_0 + a_1(x - x_i) + a_2(x - x_i)^2 + a_3(x - x_i)^3, \]

where \( a_i \) are the polynomial coefficients, expressed in terms of the second derivatives as:

\[\begin{split} \begin{aligned} a_0 &= y_i \\ a_1 &= \frac{\Delta y_i}{\Delta x_i} - \frac{\Delta x_i}{6} \left[ S''_{i + 1} + 2 S''_{i} \right] \\ a_2 &= \frac{S''_i}{2} \\ a_3 &= \frac{1}{6 \Delta x_i} \left[ S''_{i + 1} - S''_i \right] \end{aligned} \end{split}\]

Grid calculation¶

class NonuniformGridCalculator¶

Find and interpolate real numbers on a nonuniform grid.

The end points of the grid are extrapolated outward as constant values.

class UniformLogGridCalculator¶

Find and interpolate values on a uniform log energy grid.

Note that linear interpolation is applied with energy points, not log-energy points.

UniformLogGridCalculator calc(grid, params.reals);
real_type y = calc(particle.energy());

class RangeCalculator¶

Find and interpolate range on a uniform log grid.

RangeCalculator calc_range(xs_grid, xs_params.reals);
real_type range = calc_range(particle);

Below the minimum tabulated energy, the range is scaled:

\[ r = r_\mathrm{min} \sqrt{\frac{E}{E_\mathrm{min}}} \]

class SplineCalculator¶

Find and interpolate cross sections on a uniform log grid with an input spline-order.

Todo:: Currently this is hard-coded to use “cross section grid data” which have energy coordinates uniform in log space. This should be expanded to handle multiple parameterizations of the energy grid (e.g., arbitrary spacing needed for the Livermore sampling) and of the value interpolation (e.g. log interpolation). It might also make sense to get rid of the “prime energy” and just use log-log interpolation instead, or do a piecewise change in the interpolation instead of storing the cross section scaled by the energy.

SplineCalculator calc_xs(xs_grid, xs_params.reals);
real_type xs = calc_xs(particle);

class XsCalculator¶

Find and interpolate scaled cross sections.

This cross section calculator uses the same representation and interpolation as Geant4’s physics tables for EM physics:

The energy grid is uniformly spaced in \(\log E\),
Values greater than or equal to an index \(i'\) are scaled by E and are stored on a separate energy grid also uniformly spaced in \(\log E\) but not necessarily with the same spacing,
Linear interpolation between energy points is used to calculate the final value, and
If the energy is at or above the \(i'\) index, the final result is scaled by \(1/E\).

This scaling and interpolation exactly reproduces functions \( f(E) \sim a E + b \) below the \(E'\) threshold and \( f(E) \sim \frac{a'}{E} + b' \) above that threshold.

Note that linear interpolation is applied with energy points, not log-energy points.

XsCalculator calc_xs(grid, params.reals);
real_type xs = calc_xs(particle.energy());

Celeritas additionally supports two-dimensional interpolation, needed for some sampling routines.

class TwodGridCalculator¶

Find and do bilinear interpolation on a nonuniform 2D grid of reals.

Values should be node-centered, at the intersection of the two grids.

TwodGridCalculator calc(grid, params.reals);
real_type interpolated = calc({energy.value(), exit});
// Or if the incident energy is reused...
auto calc2 = calc(energy.value());
interpolated = calc2(exit);

class TwodSubgridCalculator¶

Do bilinear interpolation on a 2D grid with the x value preselected.

This is usually not called directly but rather given as the return result of the TwodGridCalculator.

Inverse sampling¶

template<class G, class C> class InverseCdfFinder¶

Given a sampled CDF value, find the corresponding grid value.

Both the input grid and the CDF must be monotonically increasing. The sampled CDF value must be in range.

Template Parameters:

G – Grid, e.g. UniformGrid or NonuniformGrid
C – Calculate the CDF at a given grid index

Construction¶

class SplineDerivCalculator¶

Calculate the second derivatives of a cubic spline.

Determining the polynomial coefficients \( a_0, a_1, a_2 \) and \( a_3 \) of a cubic spline \( S(x) \) requires solving a tridiagonal, linear system of equations for the second derivatives. For \( n \) points \( (x_i, y_i) \) and \( n \) unknowns \( S''_i \) there are \( n - 2 \) equations of the form

\[ h_{i - 1} S''_{i - 1} + 2 (h_{i - 1} + h_i) S''i + h_i S''_{i + 1} = 6 r_i, \]

where \( r_i = \frac{\Delta y_i}{h_i} - \frac{\Delta y_{i - 1}}{h_{i - 1}} \) and \( h_i = \Delta x_i = x_{i + 1} - x_i \).

Specifying the boundary conditions gives the remaining two equations. Natural boundary conditions set \( S''_0 = S''_{n - 1} = 0 \), which leads to the following initial and final equations:

\[\begin{split} \begin{aligned} 2 (h_0 + h_1) S''_1 + h_1 S''_2 &= 6 r_1 \\ h_{n - 3} S''_{n - 3} + 2 (h_{n - 3} + h_{n - 2}) S''_{n - 2} &= 6 r_{n - 2}. \end{aligned} \end{split}\]

The points \( x_0, x_1, \dots , x_{n - 1} \) where the spline changes from one cubic to the next are called knots. “Not-a-knot” boundary conditions require the third derivative \( S'''_i \) to be continuous across the first and final interior knots, \( x_1 \) and \( x_{n - 2} \) (the name refers to the polynomials on the interval \( (x_0, x_1) \) and \( (x_1, x_2) \) being the same cubic, so \( x_1 \) is “not a knot”). This constraint gives the final two equations:

\[\begin{split} \begin{aligned} \frac{(h_0 + 2 h_1)(h_0 + h_1)}{h_1} S''_1 + \frac{h_1^2 - h_0^2}{h_1} S''_2 &= 6 r_1 \\ \frac{h_{n - 3}^2 - h_{n - 2}^2}{h_{n - 3}} S''_{n - 3} + \frac{(h_{n - 3} + h_{n - 2})(2 h_{n - 3} + h_{n - 2})}{h_{n - 3}} S''_{n - 2} &= 6 r_{n - 2} \end{aligned} \end{split}\]

Once the system of equations has been solved for the second derivatives, the derivatives \( S''_0 \) and \( S''_{n - 1} \) can be recovered:

\[\begin{split} \begin{aligned} S''_0 &= \frac{(h_0 + h_1) S''_1 - h_0 S''_2}{h_1} \\ S''_{n - 1} &= \frac{(h_{n - 3} + h_{n - 2}) S''_{n - 2} - h_{n - 2} S''_{n - 3}}{h_{n - 3}} \end{aligned} \end{split}\]