.. ====================================================================== .. © (or copyright) 2019-2026. Triad National Security, LLC. All rights .. reserved. This program was produced under U.S. Government contract .. 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is .. operated by Triad National Security, LLC for the U.S. Department of .. Energy/National Nuclear Security Administration. All rights in the .. program are reserved by Triad National Security, LLC, and the .. U.S. Department of Energy/National Nuclear Security .. Administration. The Government is granted for itself and others acting .. on its behalf a nonexclusive, paid-up, irrevocable worldwide license .. in this material to reproduce, prepare derivative works, distribute .. copies to the public, perform publicly and display publicly, and to .. permit others to do so. .. ====================================================================== .. _interpolation: Gridding for Interpolation =========================== nSpiner performs interpolation on Cartesian-product grids. There are two lower-level objects: * ``RegularGrid1D`` * ``PiecewiseGrid1D`` * ``NonUniformGrid1D`` * ``FastNonUniformGrid1D`` These objects contain the metadata required for interpolation operations and have a few useful userspace functions, which are described here. Like ``DataBox``, these grid objects are templated on underlying data type, the default type being a ``Real`` as provided by ``ports-of-call``. You may wish to specialize to a specific type with a type alias such as: .. code-block:: cpp using RegularGrid1D = Spiner::RegularGrid1D; using PiecewiseGrid1D = Spiner::PiecewiseGrid1D; using NonUniformGrid1D = Spiner::NonUniformGrid1D; using FastNonUniformGrid1D = Spiner::FastNonUniformGrid1D; .. note:: In the function signature below we refer to ``T`` and ``Real`` as the underlying arithmetic data type. When constructing a ``DataBox``, you may wish to specify which interpolation object you are using. It is a template parameter. .. _grid-weight-derivatives: Interpolation weights and derivatives ------------------------------------- All four built-in grid types provide ``weightsWithGrad`` for computing the two linear interpolation weights and their coordinate derivatives in one cell lookup: .. code-block:: cpp int ix; Spiner::weights_t w, dw; grid.weightsWithGrad(x, ix, w, dw); The portable, const member function has the signature: .. code-block:: cpp void weightsWithGrad(const T &x, int &ix, Spiner::weights_t &w, Spiner::weights_t &dw) const; ``ix`` identifies the lower stencil point, and ``w[0]`` and ``w[1]`` are exactly the weights returned by ``weights(x, ix, w)``. ``dw[0]`` and ``dw[1]`` are their derivatives with respect to ``x``. For a selected interval of width ``h``, these are ``-1/h`` and ``1/h`` respectively. Custom linear grid types can implement this same interface to support :ref:`DataBox interpolation with gradients `. ``RegularGrid1D`` ------------------ We begin by discussing ``RegularGrid1D``, as the ``PiecewiseGrid1D`` object is built on top of it. Constructing a ``RegularGrid1D`` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A ``RegularGrid1D`` requires three values to specify an interpolation grid: the minimum value of the independent variable, the maximum value of the independent variable, and the number of points on the grid. These are passed into the constructor: .. cpp:function:: RegularGrid1D::RegularGrid1D(T min, T max, size_t N); Default constructors and copy constructors are also provided. Mapping an index to a real number and vice-versa ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The function .. cpp:function:: T RegularGrid1D::x(const int i) const; returns a "physical" position on the grid given an index ``i``. The function .. cpp:function:: int RegularGrid1D::index(const T x) const; returns the index on the grid of a "physical" value ``x``. The function .. cpp:function:: T RegularGrid1D::min() const; returns the minimum value on the independent variable grid. The function .. cpp:function:: T RegularGrid1D::max() const; returns the maximum value on the independent variable grid. The function .. cpp:function:: int RegularGrid1D::nPoints() const; returns the number of points in the independent variable grid. ``NonUniformGrid1D`` ---------------------- ``NonUniformGrid1D`` stores one explicitly supplied coordinate for each grid point. It is appropriate when spacing cannot be represented efficiently as regular or piecewise-regular intervals. Each grid has independent coordinate storage, so a ``DataBox`` can use a different non-uniform coordinate sequence for every axis. Constructing a ``NonUniformGrid1D`` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Constructing from a ``std::vector`` or initializer list makes an owning host copy of the points: .. code-block:: cpp NonUniformGrid1D grid({-1.0, -0.5, 0.25, 2.0}); Construction always allocates host coordinates. Use ``getOnDevice()`` when a device copy is needed: .. code-block:: cpp NonUniformGrid1D host_grid({-1.0, -0.5, 0.25, 2.0}); NonUniformGrid1D device_grid = host_grid.getOnDevice(); The coordinates must contain at least two finite, strictly increasing values. The pointer constructor borrows caller-owned memory instead and validates those requirements in the supplied execution space (the host by default): .. code-block:: cpp std::vector points = {-1.0, -0.5, 0.25, 2.0}; NonUniformGrid1D view(points.data(), points.size()); For device-resident points, pass the execution space that can access them to the pointer constructor. The supplied execution space must match the memory space of the borrowed points; using an inaccessible pointer can cause a segmentation fault. Borrowed points must remain alive and unchanged for the grid's lifetime. Like ``DataBox``, ordinary grid copies are shallow reference-style copies that also copy the ownership status, so the copy and the original both report that they own the coordinates. Finalize an owned grid exactly once, through only one of them. Every grid type also provides two explicit copy methods, which make the ownership of the result unambiguous: .. cpp:function:: void NonUniformGrid1D::shallowCopy(const NonUniformGrid1D& other) is the same as ordinary assignment, except that the result is marked ``Unmanaged``. It is a non-owning handle to the coordinates of ``other``: finalizing it never frees them, ``other`` remains responsible for them, and ``other`` must outlive it. ``shallowCopy`` may be called on host or device. .. cpp:function:: void NonUniformGrid1D::deepCopy(const NonUniformGrid1D& other) allocates independent host-owned coordinate storage and copies the coordinates from a host-resident ``other`` grid. It must not be used with a device-resident source; use ``getOnDevice()`` to make a deep device copy. Both methods require that the grid being copied into does not own memory, so that nothing is leaked. For ``RegularGrid1D`` and ``PiecewiseGrid1D``, which do not own non-trivial dynamic memory, ``shallowCopy`` and ``deepCopy`` are both trivial metadata copies. For example, an owning handle and a non-owning handle to the same coordinates can be made with: .. code-block:: cpp NonUniformGrid1D owner(points); // AllocatedHost NonUniformGrid1D view; view.shallowCopy(owner); // Unmanaged // ... view.finalize(); // no-op owner.finalize(); // frees the coordinates Mapping and interpolation ^^^^^^^^^^^^^^^^^^^^^^^^^ ``x(i)`` returns the stored coordinate. ``index(x)`` finds its lower bracketing point with a portable binary search, and ``weights(x, ix, w)`` uses the spacing of that local interval. Lookup is ``O(log N)``. Values below or above the coordinate range use the first or last interval respectively, just as ``RegularGrid1D`` does, so interpolation extrapolates linearly. ``FastNonUniformGrid1D`` ------------------------ ``FastNonUniformGrid1D`` owns the same explicit physical coordinates as ``NonUniformGrid1D`` and can add an integer lookup table that makes interval lookup ``O(1)``. The table is uniform after transforming coordinates with the first-order Ports-of-Call NQT ``asinh`` function. This provides signed logarithmic spacing far from zero and linear spacing near zero. The ``Settings`` struct supplies construction settings. A positive ``scale`` sets the transition length in the physical coordinate's units. By default it is negative, which infers the scale as the smallest coordinate magnitude at least ``PortsOfCall::Robust::SMALL()``. The resolved positive scale is available through ``scale()`` and is stored in HDF5: .. code-block:: cpp FastNonUniformGrid1D grid(points); By default, the lookup table may contain at most 32 times as many entries as the physical grid has points. If an exact table would exceed that limit, the grid transparently uses the wrapped binary search. Pass a ``Settings`` object to select a scale, policy, or different limit: .. code-block:: cpp using Settings = FastNonUniformGrid1D::Settings; using Policy = FastNonUniformGrid1D::Policy; FastNonUniformGrid1D grid( points, Settings{ .scale = scale, .policy = Policy::RequireFast, .max_lookup_ratio = 16}); ``Automatic`` permits fallback, ``RequireFast`` fails construction if the table cannot fit, and ``ForceBinary`` skips the table. The selected mode can be inspected with ``usesFastLookup()`` and ``lookupSize()``. A host-owned grid can be reconfigured later, including after HDF5 loading: .. code-block:: cpp auto settings = grid.settings(); settings.scale = -1.0; // infer again from the physical coordinates settings.policy = FastNonUniformGrid1D::Policy::ForceBinary; grid.reconfigureLookup(settings); Reconfiguration follows the same explicit ownership convention as ``finalize()``: do not reconfigure an owner while shallow aliases depend on its lookup table. Reconfigure host storage before making device copies. Ordinary copies remain shallow. ``shallowCopy()`` makes a non-owning handle to both the coordinates and the lookup table. ``deepCopy()``, ``getOnDevice()``, binary serialization, and HDF5 otherwise follow the ``NonUniformGrid1D`` lifecycle. HDF5 stores the physical coordinates and resolved construction settings, then rebuilds the derived lookup table when loading. Coordinate access is read-only because changing a coordinate would invalidate the table. The ``PiecewiseGrid1D`` ------------------------ A ``PiecewiseGrid1D`` is a non-intersecting, contiguous, ordered collection ``RegularGrid1D`` s. It can be used to construct grids with non-uniform spacing, so long as the grid spacing is piecewise constant. The maximum number of ``RegularGrid1D``s that can be used to construct a ``PiecewiseGrid1D`` is a compile-time parameter (default is 5). You can specify a different value with, e.g., .. code-block:: cpp // Maximum number of "pieces" in a grid = 10 using PiecewiseGrid1D = Spiner::PiecewiseGrid1D; Constructing a ``PiecewiseGrid1D`` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A ``PiecewiseGrid1D`` is constructed from either a ``std::vector`` or a ``std::initializer_list`` of ``RegularGrid1D`` s. For example: .. code-block:: cpp // Initialize the regular grids // Note that the start and end points match // for each consecutive pair of grids. // g1 ends when g2 starts, etc. Spiner::RegularGrid1D g1(0, 0.25, 3); Spiner::RegularGrid1D g2(0.25, 0.75, 11); Spiner::RegularGrid1D g3(0.75, 1, 7); // Build the piecewise grid. The double bracket notation // is an "initalizer list" and is very convenient, // as it is a C++ language feature. Spiner::PiecewiseGrid1D h = {{g1, g2, g3}}; Default constructors and copy constructors are also provided. Index Mapping with ``PiecewiseGrid1D`` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A ``PiecewiseGrid1D`` has all the same functionality as ``RegularGrid1D``, but it automatically uses the relevant underlying grid spacing. The function .. cpp:function:: T PiecewiseGrid1D::x(const int i) const; returns a "physical" position on the grid given an index ``i``. The function .. cpp:function:: int PiecewiseGrid1D::index(const T x) const; returns the index on the grid of a "physical" value ``x``. The function .. cpp:function:: T PiecewiseGrid1D::min() const; returns the minimum value on the independent variable grid. The function .. cpp:function:: T PiecewiseGrid1D::max() const; returns the maximum value on the independent variable grid. The function .. cpp:function:: int PiecewiseGrid1D::nPoints() const; returns the number of points in the independent variable grid. Developer functionality ------------------------ All grid types implement the resource-management interface used by ``DataBox``: .. cpp:function:: std::size_t Grid::serializedSizeInBytes() const; .. cpp:function:: std::size_t Grid::dynamicMemorySizeInBytes() const; .. cpp:function:: std::size_t Grid::dumpDynamicMemory(char *dst) const; .. cpp:function:: std::size_t Grid::serialize(char *dst) const; .. cpp:function:: std::size_t Grid::deSerialize(char *src); .. cpp:function:: std::size_t Grid::setPointer(char *src); .. cpp:function:: void Grid::shallowCopy(const Grid& other); .. cpp:function:: void Grid::deepCopy(const Grid& other); .. cpp:function:: Grid Grid::getOnDevice() const; .. cpp:function:: void Grid::finalize(); These operations are currently trivial for ``RegularGrid1D`` because it contains only inline data. ``PiecewiseGrid1D`` applies them recursively to its component grids. ``NonUniformGrid1D`` uses them to manage its coordinate array. This interface allows grid types to own dynamically allocated host or device data without requiring grid-specific resource handling in ``DataBox``. ``DataBox::finalize`` does not call ``Grid::finalize``. ``DataBox::finalizeGrids`` does, for every grid in the ``DataBox``. ``DataBox`` copies its grids with ordinary assignment when it makes a new handle to the same ``DataBox``. When it builds a different ``DataBox`` from another's metadata, it uses ``shallowCopy`` or, for deep copies, ``deepCopy``. See :ref:`Semantics and Memory Management` for details. ``serialize`` writes the inline grid object followed by its dynamic memory. ``dumpDynamicMemory`` writes only the latter, allowing a grid embedded in a ``DataBox`` to avoid serializing its inline bytes twice. The binary serialization methods have the same transient, build-dependent compatibility limitations as ``DataBox`` serialization. See :ref:`the DataBox serialization documentation ` for details. Generative AI was used to assist with modifications to this page.